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Abstract 

I review the current state of numerical simulations of stellar feedback in the context of star formation at scales ranging 
from the formation of individual stars to models of galaxy formation including cosmic reionisation. I survey the wealth 
of algorithms developed recently to solve the radiative transfer problem and to simulate stellar winds, supernovae and 
protostellar jets. I discuss the results of these simulations with regard to star formation in molecular clouds, the 
interaction of different feedback mechanisms with each other and with magnetic fields, and in the wider context of 
galactic- and cosmological-scale simulations. 
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1 Introduction and scope of this review 


The formation of stars is arguably the most important process in astrophysics, impacting virtually every theoretical 
and observational subfield. Despite its prominence and decades of intensive study, there is still much about star 
formation that is not understood. One reason for this slow progress is the fact that the conversion of gas to stars is a 
non-linear process. There are a variety of reasons for this, such as the non-linearity of the self-gravitational forces 
which lead to the collapse of individual stars, but the issue of most interest here is that the rate of star formation is 
influenced by stars themselves via feedback. 

There are many other phenomena that affect the rate and morphology of star formation, such as mergers 
between galaxies (e.g. flWhitmore et al., 2010| , [Ho pkins et al., 20TTj ), collisions between molecular clouds 
(e.g. jjFurukawa et al., 2009| , QTasker and Tan, 20091 ) or the passage of spiral shocks through galactic discs (e.g. 
IjMeidt et al., 20131 , jjBonnell et al., 20131 ). However, while a given star-forming region may or may not have 
experienced these particular perturbations, feedback from the stars themselves is, by definition, always present, and it 
is this broad range of processes that are the focus of this review 

Stellar feedback has been invoked, with varying degrees of success, to solve a wide range of issues and problems 
in astrophysics, including the slow and inefficient star formation observed in molecular clouds and on galactic and 
cosmological scales, the triggering of star formation, the formation of disc galaxies and the suppression of excess 
dwarf galaxy formation in cosmological simulations. A glance at almost any image from HST, Spitzer, Herschel, 
WISE and many ground-based images reveals that the structure of the ISM is riddled with bubbles, shells, pillars and 
outflows, none of which can be explained without invoking feedback. Since one of the main purposes of astrophysical 
simulations is to help explain what is observed in the Universe, it is clear that feedback is a critical component of such 
simulations, and of any general model of star formation. 

More specifically, this is a review of numerical simulations of feedback. A self-gravitating pure-hydrodynamics 
problem would already by of sufficient complexity to require the use of high-dimensional computer simulations. 
The inclusion of additional physical processes, particularly the transfer of radiation, only makes the problem more 
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complex and the necessity of using simulations all the greater. However, progress in this field is very rapid and 
the existence of a varied set of mutually-interacting feedback processes has resulted in a bewildering number of 
recent studies. A timely summary will therefore be of beneht to specialists and non-specialists alike. In this 
review, the fundamental physical processes will be only briefly rehearsed; a detailed overview can be found in 
OKrumholz et al., 2014| . The spotlight will instead fall on the algorithms that have been written to model them, 
and of simulations which have been performed including them, and what we have learned from these simulations. 
Numerical studies of feedback have a long and rich history (e.g. | Tenorio-Tagle, 19791, | Tenorio-Tagle et al., 19851, 
[Yorke et al., 1989) , |Garcia-Segura and Franco, 1996| ). No single review could encompass all this work, and this 
article will concentrate on articles published for the most part in the last ten years, and on two- or three-dimensional 
simulations. 

Although feedback in the galactic context will be discussed, this review deals with stellar feedback and, despite its 
evident importance, AGN feedback will not be covered. Interested readers are referred to QFabian, 2012| . Similarly, 
since we are here concerned with the connection between stellar feedback and star formation, the focus will be on 
feedback from low- and intermediate-mass proto- and pre-main stars, and on O-stars, whose entire life-cycle is 
comparable to the lifetimes of GMCs. Planetary nebulae and feedback from accretion onto compact objects will not 
be examined, and readers are instead directed to QBalick and Frank, 2002| |. For the most part, this review also does 
not cover feedback from Population III stars, an area of research which has grown substantially in recent years, and 
which is eloquently reviewed by [Greif, 2015| . 

Magnetic fields are not commonly regarded as a type of feedback and will not receive dedicated attention here. 
However, the presence of a magnetic field will likely alter the response of a fluid to some or all of the feedback 
mechanisms under consideration and several authors have performed simulations including both feedback and 
magnetic fields. This work will be discussed, but the algorithms used to model the magnetic held will not be 
described. 

The structure of this article is as follows. Section 2 gives a brief introduction to the major feedback mechanisms, 
namely photoionisation, stellar winds, supernovae, accretion heating, radiation pressure and protostellar jets. Section 
3 briefly introduces the major classes of astrophysical hydrodynamics codes - particle-based schemes such as SPH, 
grid-based schemes such as AMR, and the new generation of moving-mech codes. Section 4 surveys the algorithms 
used for modelling radiation transport, winds, supernova and jets. Section 5 discusses the science which has been done 
with the codes described in Section 4 with reference to particular astrophysical problems, including the fragmentation 
and destruction of molecular clouds and the formation and evolution of spiral and dwarf galaxies. Section 6 contains 
a short summary and outlook for the future. 


2 Brief introduction to stellar feedback physics 

Stellar feedback involves the insertion of matter, momentum and energy from stars into the surrounding fluid, from 
which the stars may also still be accreting gas. In terms of material, momentum and energy emitted per star, massive 
OB-type stars far outweigh their lower-mass brethren in importance. In clouds where there are no O-stars (either 
because the cloud mass is too small to support massive star formation, or because there has not been time for O-stars 
to form), feedback from low- and intermediate-mass stars in the form of jets and outflows and the conversion 
of gravitational potential energy to heat are dominant processes. On larger lengthscales and longer timescales, 
encompassing the formation and evolution of galaxies and cosmic star formation, it is again the O-stars that dominate 
the stellar contribution to feedback, but the other smaller-scale processes may still have influence, since they help 
determine the environments in which the O-stars are born. 


2.1 Photoionisation 


The physics of photoionisation were first elucidated in detail by |Str6mgren, 19391 and an excellent modern intro¬ 
duction can be found in [Osterbrock and Ferland, 2006) . If a source of Qh ionising photons per second ignites in a 
cloud with initial number density no atoms cm~^, the number density of ions ni will also initially be equal to no. If 
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the cloud is pure hydrogen and overall electrically neutral, the electron number density Ug must equal the ion number 
density, so we have rig = rii = no. 

As the ionisation front moves outwards, more photons are consumed by recombinations behind it. The recombi¬ 
nation rate per unit volume is given by aUgUi = auQ. Recombinations directly to the ground state produce photons 
which are able to ionise another atom elsewhere in the nebula. However, since a recombination of this kind quickly 
produces another ion, their overall effect can be neglected to first order, which is known as the ‘on-the-spot’ (OTS) 
approximation. Only recombinations to states other than the ground state consume stellar photons, and this rate is 
given by aBn^, ctb being the recombination rate to all states above the ground state. 

Eventually, the total recombination rate behind the ionisation front equals the rate at which the source produces 
photons, so no more neutral gas can be ionised. This state is known as the Stromgren sphere, described by the 
Stromgren radius Rg, satisfying 


Rs = 


3Qh \ ^ 

dTrasriQ J 


( 1 ) 


The ionised gas inside the Stromgren sphere has a temperature determined by the equilibrium of heating and cooling 
processes. The main heating process is the absorption of stellar photons. In a nebula with solar metallicity, the 
main coolants are optical line emission from metals, with some contribution from free-free emission. These process 
equilibrate at 8 000 - 10 OOOK, although this temperature may be higher in low-metallicity environments where line 
cooling is less efficient. Since the temperatures of neutral clouds lie in the range 10-100 K, the Stromgren sphere will 
be vastly overpressured and will expand supersonically. This process was first studied by [Spitzer, 1978[ , who derived 
the well-known relation for the radial evolution of the shock driven by the expanding Hll region 


R{t) = Rs 



( 2 ) 


Later, [Franco et al., 1990) considered the case of HII regions expanding in clouds with radial density gradients 
p{r) oc r“". They showed, in particular, that if n > 3/2, the ionisation front will inevitably overtake the shock front 
and ionise the whole cloud, regardless of its mass. Early 2D numerical work by [Garcia-Segura and Franco, 1996| 
found that expansion of ionisation fronts in density gradients, and in uniform clouds, was accompanied by the forma¬ 
tion of fingers of dense neutral material reaching into the HII region, reminiscent of observed pillars/elephants’ trunks. 
They interpreted these structures in terms of the generic ionisation front instability analysed by [Giuliani, 1979|. 


2.2 Main-sequence line-driven winds 

The powerful fluxes of energetic photons emitted by OB stars are able to accelerate line-driven winds in their 
atmospheres. Material leaves the surface of the star at a velocity which asymptotically approaches the wind terminal 
velocity Voo, of order 10^ km s“^ for a main-sequence O-type star ( [Earners and Cassinelli, 1999| ). The wind mass 
fluxes M for such stars are typically ^ yr“^ but can approach ^ IO^'^Mq yr“^. 

The effect of the wind depends sensitively on the thermodynamics of the shocked gas inside the wind bubble. 
The extremal assumptions are (i) that all the mechanical energy is retained by the bubble and the expansion is 
adiabatic, or pressure driven, in which case R{t) = (L/po) = , and (ii) that cooling is maximally efficient, whence 

R{t) = [2Muoo/(3po)]Note in both cases the very weak dependence on the both the stellar properties and 
the density of the background medium. The reality is much more complex and lies somewhere between these 
extremes. More sophisticated models were first computed in ID by [Castor et al., 19751 [Weaver et al., 1977[ who 
examined the influence of thermal conduction between the hot shocked wind and the cool swept-up ISM and 
the corresponding evaporation of ISM material into the wind bubble, which rapidly comes to dominate its mass. 
[Koo and McKee, 1992a| and [Koo and McKee, 1992b| discuss in exhaustive detail the realistic case of partially 
radiative bubbles. 

Recently, [Rosen et al., 2014) have attempted to evaluate the importance of these processes observationally by 
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totting up the energy inserted by winds and that lost by various physical processes (radiative cooling, mechanical 
work, thermal conduction, gas-grain interactions) in four star-forming regions - 30 Doradus, Carina, NGC 3603 and 
M17. They concluded that radiative cooling and mechanical work are unimportant, except in the case of M17 where 
work done on the cold gas can account for 38% of the injected wind energy. Adding thermal conduction and gas-dust 
cooling can account for the remainder of the input energy, but only if rather extreme assumptions are made. They 
instead infer that large fractions of the wind energy is lost via bulk leakage of hot wind gas, or small-scale mixing 
with the cold gas. 


2.3 Stellar evolution 

The lifetimes of the most massive stars are comparable to or shorter than the inferred lifetimes of GMCs, so at least 
some of these stars are likely to enter the hnal stages of their lives while embedded in their natal clouds. The effects of 
exotic evolutionary phases such as LBV and WR stars have generally not been considered in GMC-scale simulations 
owing to the expense of simulating the long timescales involved. 

The WR stage profoundly alters the properties of the wind of a massive star. The mass loss rate increases 
dramatically to ^ 10 “^ — yr“^ as the star sheds its outer layers, and the terminal velocity correspondingly 

declines to ^ 100km s“^. A thorough introduction to these and other kinds of stellar wind can be found in 
[Tamers and Cassinelli, 199^ 

Understanding the effects of photoionisation and winds clearly relies on knowing how the ionising photon 
luminosity and wind mass loss rate and terminal velocity vary as a function of stellar age and mass. There are several 
observational studies of this issue, such as [Smith, 2006[ or [Weidner and Vink, 2010) . 

Photoionisation and winds have traditionally been the most popular feedback mechanisms, perhaps since 
their effects are readily observable as ~ lOpc bubble structures in atomic emission lines, radio continuum 
and infrared dust emission. Several authors have considered which of the two should be more important 
(e.g. [Capriotti and Kozminski, 2001| [Matzner, 2002) ), generally concluding that expanding HIIRs are more damag¬ 
ing. 


2.4 Supernovae 

When massive (> 8 M 0 ) stars exhaust their core hydrogen, a chain of events ensues which eventually results in a 
supernova explosion. The timescale on which hydrogen exhaustion occurs depends on the stellar mass. For a IOMq 
star, it is «30 Myr, but for the most massive stars, it may be as short as Ri 5Myr, comparable to or shorter than the 
lifetimes of GMCs. The supernova explosion results in the ejection of ^10 Mq of metal-enriched material at speeds 
of Ri 3 X lO^km s“^, carrying approximately lO^^erg of total energy. In the classic problem of the point deposition 
of energy in a uniform medium, the supernova remnant passes though a brief ‘free-expansion’ phase until the mass 
swept up becomes comparable to the ejecta mass, before entering the adiabatic Sedov-Taylor phase, during which the 
blast wave radius evolves with time as 


R{t)=^(- 

\Po 




(3) 


where (3 is an order-unity constant. The Sedov-Taylor phase ends when the cooling timescale becomes shorter than 
the expansion timescale and the supernovae remnant enters the radiative phase. In reality, the evolution of a supernova 
blastwave is likely to be much more complex, since it in reality encounters the interior of a wind bubble/HII region, 
rather than a smooth ambient medium. A comprehensive review of the expansion of astrophysical shock- and blast- 
waves, including the cases of nonuniform background density helds, is given by [Ostriker and McKee, 1988| . 

Much of our understanding of the evolution of massive stars is predicated on the assumption of single stars 
evolving in isolation, owing to the extreme difficulty of assuming otherwise. However, a large fraction (45-75%) of 
massive stars are members of binaries (usually with other massive stars), most of which are sufficiently close that the 
evolution of the two members is expected to be influenced by mass exchange or outright mergers (see [Langer, 2012) 
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for a recent review of massive single and binary star evolution). How binary evolution would affect the production of 
ionising photons, the emission of winds, and the timing and energy of supernovae is a field with a lot of work still 
ahead of it. 


2.5 Accretion heating 

In the early stages of formation, all protostars are accreting gas from their host GMCs via an envelope and a disc. 
As the gas falls into the potential well of the protostar and is eventually accreted by it, gravitational potential energy 
is converted to heat, which the protostar radiates away. Although not readily absorbed by the gas, this radiation is 
absorbed by dust. If the gas density is large enough, thermal coupling of the gas and dust may then transfer heat to 
the gas, effectively coupling the radiation held to the gas. QOffner et ah, 20091 point out that accretion luminosity 
dominates the energy budget of GMCs unless and until O-stars are born. 


2.6 Radiation pressure 


All photons emitted by all stellar objects transfer momentum, as well as thermal energy, to the surrounding gas and 
dust. The accretion heating process described above also exerts radiation pressure forces on dust grains, although 
these are likely to be small and less important than the thermal pressures generated by dust and gas heating. However, 
radiation pressure from much more luminous massive stars is likely to be dynamically important in very dense clouds 
(e.g. [Krumholz and Matzner, 2009][Murray et ah, 2010|[Fall et ah, 2010| |). 

Although a conceptually simple idea, radiation pressure is difficult to compute in practice. The momentum carried 
by the photons of a source of luminosity Lis L/c, giving (for a point source) a radiative momentum flux at radius r of 
L/( 47 rr^c). However, to compute the radiation pressure, one needs to know how much of this momentum is actually 
absorbed by the gas. Clearly, if the gas is optically thin, the absorbed momentum can be zero, even if the radiative 
momentum flux is very large. [Krumholz and Matzner, 2009] introduce a parameter /trap which encapsulates this 
uncertainty, so that 


Pra,d{r) 


/trap 


dTTr^C 


(4) 


If /trap = 1, all the emitted photons are absorbed once before escaping. However, if the shell is moderately optically 
thick, photons are likely to interact several times, depositing more momentum, before escaping the shell and /trap 
exceeds unity. Clearly, computing /trap self-consistently is a very difficult radiative transfer problem. 


2.7 Jets 

As well as thermal feedback, accretion also drives emission of stellar jets - collimated high-velocity outflows 
emerging bidirectionally along the stellar rotation axis (for a recent view, see [Frank et al., 20141 ). The origin of jets 
is likely a magneto-hydrodynamic interaction between the star and its accretion disc, producing a magnetic field 
configuration which acts like a particle accelerator or collimator. The details of this process are still much debated 
and I will not address them here. From the point of view of view of feedback, it is useful to know that initial jet 
velocities fall in the range 100-1 000 kms“^, while material swept up by the jet bowshocks has typical velocities 
in the range 1-30 kms“^. The mass-loss rate via jets is typically ~lO“®M 0 yr~^. The shocks produced by jets are 
highly radiative, so to a good approximation they may be considered as sources of momentum only. 


3 Brief introduction to astrophysical fluid dynamics codes 

Star formation takes place in the interstellar medium (ISM), the thin and usually hot gas which occupies much 
of the volume of most galaxies. The mean free paths of ions, atoms and molecules in the ISM tend to be small 
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compared with the sizes of the structures which they belong to. It is therefore reasonable to approximate the ISM 
as a smoothly-varying fluid. However, in order to model the behaviour of a fluid on a computer, it is necessary to 
discretise it in some way into individual elements. How the discretisation is done inevitably has some effect on 
how physical processes are treated. A very brief summary of the three main types of astrophysical hydrodynamics 
codes and their advantages and disadvantages is therefore necessary before specific feedback algorithms are discussed. 
Since this review covers simulations of star formation, particular attention will be paid for each type of code to the 
way in which the formation of individual stars is modelled. 


3.1 Grid codes 


Grid codes break fluids up into volume elements which All the space inside a set of boundaries delimiting the 
computational domain. The simulation is evolved by calculating the forces on the fluid in each volume element 
and moving material to or from that element’s adjoining neighbours. Matter which crosses one of the boundaries is 
either destroyed (open or outflow boundaries), sent back into the domain as though bouncing off a wall (reflective 
boundaries), or reinserted at the opposite side of the domain (periodic boundaries). Material may also be created at the 
boundaries and allowed to flow into the domain (inflow boundaries). Some codes make do with a single fixed grid, but 
most allow a given grid cell to be subdivided or refined if higher resolution is required at that location, or incorporated 
into larger grid cells if the resolution at that place is higher than needed. Codes such as FLASH [ Fryxell et ah, 2000) 
and RAMSES which are able to do this on the fly are known as Adaptive Mesh Refinement codes. 

The hydrodynamic equations themselves are solved by a wide variety of methods. Finite difference methods 
discretise the differential equations connecting quantities in adjacent cells. Finite volume methods integrate quantities 
over the volumes of the grid cells to compute fluxes between them, often by solving the Riemann problem. Exhaustive 
reviews of these various methods can be found in many textbooks, e.g. [Bodenheimer et ah, 2007| . 

The advantages of grid codes include being able to use virtually any criterion to decide when and where to refine or 
dereflne the grid, and that the mass contained in individual grid cells may become very small or very large, so that very 
large dynamic ranges in the masses of objects are possible. Disadvantages are that some form of boundary condition 
must be specified, which limits the volume that can be studied, all grid cells must contain non-zero quantities of gas, 
so that some computational power may be wasted on simulating regions where very little is happening, and that fluid 
advection is almost inevitably slightly more efficient along the principal grid axes, leading to artefacts (often known 
as ‘carbuncles’). Tracing fluid flows in grid codes can be done by advecting passive scalars, but this is somewhat 
cumbersome and cannot be used to trace completely arbitrary flows. Modelling gravity in grid codes is non-trivial 
and is usually done by solving Poisson’s equation using multigrid methods, although there are also implementations 
of tree algorithms similar to those used in particle-based codes. In addition, grid codes are not Galilean invariant and 
moving objects through the grid at speeds in excess of the local sound speed can lead to problems such as unphysical 
diffusion. 

If it occurs that the density of any particular grid cell becomes so high that its integration time becomes 
prohibitively short, some of the mass in the cell may be converted into a Lagrangian sink particle which is allowed to 
move between grid cells and to accrete further material. [Federrath et ah, 2010| describe in detail their implementation 
of sink particles in FLASH. As with SPH codes, the first criterion is that the density of the grid cell in question must 
exceed a given threshold. Six tests are then applied for all cells within a sink accretion radius of the dense cell: (i) the 
cells must all be on the highest allowed AMR refinement level; (ii) the flow at that location must be converging; (iii) 
the densest cell must lie at a local potential minimum; (iv) the gas within the accretion radius must be Jeans unstable; 
(v) the gas within the accretion radius must be bound; (vi) the volume must not overlap the accretion volume of a 
pre-existing sink. Once a sink is created, it may later accrete gas above the threshold density in the grid cell in which 
it finds itself, provided that gas is bound to it. 


3.2 Particle codes 

Particle codes discretise fluids into mass elements with a total mass equal to that of the whole fluid. The fluid is 
evolved by calculating the forces on each particle and moving the particles relative to one another. Since the particles 
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can in principle go anywhere, it is not obligatory to have boundaries of any kind in a particle simulation, although 
reflective and periodic boundaries are commonly used. Open boundaries at which particles are destroyed, or inflow 
boundaries where new particles are inserted are also possible but uncommon. 

By far the most popular particle method is Smoothed Particle Hydrodynamics (SPH, e.g. | Springel, 2010b |). SPH 
codes represent fluids as particles whose masses are smoothed or smeared out over a volume called the smoothing 
kernel. The kernel is spherical, but the smoothing is done such that most of the mass is concentrated near the particle 
centre. The radius of the smoothing kernel is continuously recomputed so that it contains the centres of, on average 
(or, in some codes, exactly) a given number (usually Ri 50) of ‘neighbour’ particles. This is achieved in many codes 
(e.g. SEREN, QHubber et al., 20lTl ) by iteratively solving for each particle i the expression ^constant (i.e. the 
fixed particle mass), as described in[ Springel, 2010b) . Fluid quantities such as density are computed as averages over 
a particle and all of its neighbours. The requirement that the number of neighbours be fixed automatically results in 
the smoothing kernels, which set the resolution of an SPH code, being smallest where the gas density is highest. 

Gravitational forces could in principle be computed directly between particles, but this scales extremely poorly 
as the number of particles increases so is rarely if ever used. Many codes reduce the expense of computing the 
gravitational forces by grouping particles into a tree structure, and computing gravitational forces using tree nodes 
instead of individual particles, provided that the tree node subtends a sufficiently small angle at the location where the 
forces are to be computed. Alternatively some codes use the particle-mesh method where the particle densities are 
converted to densities on a mesh and gravitational forces are computed by solving Poisson’s equation. 

Not needing boundaries, the relative ease of computing self-gravitational forces, and the possibility of having 
parts of the computational domain genuinely empty are the main advantages of particle codes. Since each particle 
has a unique and preserved identity, it is also trivial to trace fluid flows in particle simulations, which is very helpful 
when, for example, studying triggered star formation or pollution by supernovae. Particle codes generally do not treat 
shocks as well as grid codes, and they are restrictive in the sense that only the mass density can be used to control the 
resolution and other quantities cannot be refined on as needed. Traditionally, particle code also have problems dealing 
with contact discontinuities ([ Agertz et al., 2007) ), although there are now several solutions to this issue available 
(e.g. ||Read et al., 2010||Saitoh and Makino, 20131 ). 

Since they are Lagrangian codes, implementing sink particles in SPH schemes is somewhat more straightforward 
than in grid codes, since the sinks can be treated like gas particles except that they do not feel or exert pressure forces. 
An early implementation was described by [Bate et al., 1995| . A density threshold is defined and any gas particles 
exceeding this threshold, along with their neighbours, are considered for sink formation. Four criteria must be met; 
(i) the ratio a of thermal to gravitational potential energy of the group must be < 0.5 (ii) the sum of a and fl, the ratio 
of rotational energy to gravitational potential energy must be < 1 (iii) the total energy of the group must be negative 
(iv) the divergence of the acceleration must be negative. If these tests are passed, a sink is created with the total mass 
and momentum of the seed gas particles. 

Accretion onto the sink is achieved by assigning it an accretion radius and testing particles which pass within it. 
Particles which are bound to the sink with a specific angular momentum less than that required to form a circular orbit 
at the accretion radius are accreted. A much more sophisticated SPH sink particle algorithm was recently presented 
by QHubber et al., 20131 . The creation criteria are again that the gas particle being considered for promotion must 
have a density exceeding a threshold, the putative sink particle would not overlap any pre-existing sinks, it must sit 
at a local potential minimum, and the candidate particle’s density must be such that it is smaller than the Hill sphere 
defined by itself and any pre-existing sink. 

Once created, their sinks do not immediately accrete gas particles entering the accretion radius (which they term 
the ‘interaction zone’). Instead they are added to an interaction list (from which they may be struck off if they exit 
the interaction zone) and are gradually accreted over a physically-motivated timescale, while still being permitted to 
interact with other SPH particles. The smooth accretion and the continued interaction with gas particles outside the 
sink interaction zone, particularly in respect of angular momentum transfer, results in more physically-motivated and 
robust sink behaviour. 

In practice, most of these conditions are usually dropped in large-scale simulations where sink masses are too big 
for them to considered as single stars. In these cases, sink particles are often created simply from particles whose 
densities exceed a threshold. 
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Figure 1: Simple illustration of the ideas behind the major types of astrophysical hydrodynamics code. The left panel 
shows how an SPH code represents a fluid, with red dots being individual particles carrying a mass m and specific 
internal energy u. Each particle has a velocity v and feels forces from other particles. Denser regions of gas are 
represented by higher concentrations of particles. The size of a particle is defined by the radius 2h which encloses a 
constant number of neighbours (here, only 4 for clarity). The middle panel shows an AMR representation of a fluid 
with only three levels of refinement. The gas in each cell has a density p, velocity v and internal energy e. Cells 
exchange mass, momentum and energy with their neighbours. The right panel depicts a moving mesh code. The fluid 
mass is again represented by particles, shown as red dots, and the fluid volume is partitioned around these by a Voronoi 
tessellation. Cells exchange matter, momentum and energy, but the particles move with the fluid flow. 


3.3 Moving-mesh codes 


The most recent additions to the stable of astrophysical fluid codes are the moving-mesh codes (e.g. 
fGaburov and Nitadori, 20lT| Hopkins, 20151). The most widely used to date is AREPO (| Springel, 2010a[ ) which 
solves the hydrodynamical equations using a finite-volume Godunov method on a 3D Voronoi mesh dynamically 
created around a population of Lagrangian particles which follow the fluid flow. The code is then in some sense a 
hybrid between traditional grid- and particle-based codes, and shares most of the advantages and few of the draw¬ 
backs of both alternatives. Moving-mesh codes are Galilean-invariant like SPH codes but capture shocks and contact 
discontinuities as well as grid codes. They also share the ability to refine or dereflne their resolution based on arbitrary 
conditions by splitting or merging their Lagrangian tracer particles. Sink particles can be implemented in ways anal¬ 
ogous to those employed by grid-based codes, in which sinks remove mass from grid cells. The only disadvantages 
of moving-mesh codes are their relative complexity and novelty. Figure [T] illustrates schematically the differences 
between these three types of code. 


4 Feedback algorithms 

This section briefly surveys some of the algorithms used to model stellar feedback mechanisms. The focus is on the 
algorithms themselves and the assumptions that underlie them. The results gained from using them will be discussed 
in a later section. 


4.1 Radiative transfer algorithms 

During and after their formation, stars - even low-mass objects - are strong sources of radiation which deposit energy 
and momentum into the surrounding gas. The interaction of radiation and matter is an immensely difficult problem 
to solve and a great deal of effort has been expended on it. The summary given here is of necessity brief - a more 
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detailed and wide-ranging review can be found in QTrac and Gnedin, 201 1| . The two main processes of interest here 
are the radiation emitted mostly in the infrared by protostars (deriving from the loss of gravitational potential energy 
by accreting material and by the contracting protostar itself, and sometimes due to the burning of deuterium before 
hydrogen burning gets underway), and emission of ultraviolet photoionising photons by massive stars as they quickly 
settle onto the main sequence. 

In common with self-gravitational forces, radiative transfer (RT) in principle allows every part of the computa¬ 
tional domain to communicate with every other part. However, the latter issue is much worse, since the communication 
between any two fluid elements depends on all the intervening material through which any radiation must pass, which 
of course does not apply to gravitational forces. This problem is in general too demanding to be solved explicitly but 
numerous physical and numerical approximations have been made to render it tractable. 


4.1.1 Ray-tracing methods 

The most intuitive approach to RT computations is ray-tracing, that is, drawing lines from radiation sources to target 
fluid elements and solving the radiation transfer equation along them. We shall limit ourselves to the consideration of 
unpolarised radiation. The time-dependent RT equation takes the form 

oT” - €i/ (5) 

C Ot 

with 4 the specific intensity at frequency n a unit vector pointing in the direction of radiation propagation, the 
emissivity of the medium and the specific absorption coefficient. This is an equation in seven variables and solving 
it is a formidable problem. This is particularly true if large frequency ranges are of interest, for which the frequency 
dependence of the emissivity and opacity is likely to be significant and the problem has to be solved for many different 
values of v. It is common to avoid this issue by computing effective average emissivities and opacities, commonly 
referred to as the ‘grey’ approximation. 

It is also common to simplify Equation|5]by various approximations. If the time-dependence can be neglected, 
equivalent to finding a radiative equilibrium, the time-independent RT equation results: 

n.V4 = Cl/- Ki/p4- (6) 

Furthermore, it is often the case that the radiation held is dominated by a small number of very bright sources (e.g. 
stars), and that the emissivity of the gas may be neglected, yielding 

n.V4 = -K^pU. (7) 


Several authors have made use of ray-tracing algorithms to attack the problem of photoionisation, differences 
being mainly in how they choose to cast their rays. Under the OTS approximation, any number of independent rays 
may be drawn emanating from an ionising source and the thermodynamic state of the gas can be found by locating the 
ionisation front along each ray using a generalisation of Equation [T] If the radius of the ionisation front is a function 
of direction Rif{ 9, (p), one can write 


Oh 

47r 


r'=RiF{e,(j>) 


i{r, 9, (p) 


2 /2 
asr 


dr'. 
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( 8 ) 


I Kessel-Deynet and Burkert, 2000) and ||Dale et ah, 2007cl using SPH codes defined rays connecting the ionising 
source to all active gas particles (a similar method was used by [Johnson et ah, 2007) , except they first constructed 
a spherical grid with 10® rays, each divided into 500 radial segments). Particles’ neighbours are tested to And the 
one closest (in an angular sense) to the ray leading back source. This process is repeated until the source is reached, 
generating a list of particles along the ray, whose densities are then used to calculate the integral in Equation This 
can be used in a time-independent way to locate the ionisation front assuming ionisation equilibrium and heat the 
gas behind the ionisation front. However, [Dale et ah, 2007cl use it to compute the photon flux at each particle to 
determine whether it is sufficient to keep an ionised particle in that state, or to (partially or completely) ionise a neutral 
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particle, during the current timestep. 

This algorithm was modified in [Dale and Bonnell, 201 1| and JDale et ah, 2012bl to allow for multiple sources 
ionising the same Hll regions. In the former paper, this was achieved by identifying all particles illuminated by more 
than one source and dividing their recombination rates by the number of sources illuminating them. The solution for 
the radiation field is iterated until the number of ionised particles converges. In the latter paper, a more sophisticated 
approach was adopted where the total photon flux at each particle is evaluated and the fraction of the recombination 
rate that each source is expected to pay for at a given particle is set by the flux striking it from that source as a fraction 
of the total. These methods give similar results in practice. 

In the Eulerian HERACLES code by [Tremblin et al., 2012b| , solve a differential form of Equation [8] taking the 
photon flux F as the variable of interest, writing 

= Qh - nHoCr^", (9) 

dr 

where a is the absorption cross section to ionising photons and uho is the density of neutral hydrogen atoms. This 
equation is then complemented with equations describing photochemistry. The ionisation fraction x is /nn and 
riB. — + nug. The number of photons absorbed in a grid cell of volume dl4eU is given by the flux of photons 

entering the cell through the surface dA multiplied by the probability of absorption dP = crriHods, with ds being the 
pathlength along the ray. This allows one to write 


d(a;nH) = dn^ - dnn^ec: 

(10) 

dnH,ec = asx'^nlidt, 

(11) 

rn/ \ dP ^ , dAds 

(12) 


where the last term in the last equation is a geometrical dilution factor. Once these equations have been solved, the 
heating due to the absorption of ionising photons and the cooling due to recombinations can be computed. 

QGritschneder et al., 2009a| modelled the propagation of plane-parallel ionising radiation in SPH. Rather than 
drawing rays to every particle, they used an adaptive scheme, casting a small number of rays along the photon 
propagation direction, and recursively refining them into four subrays up to five times at locations where the 
separation of the rays exceeded the particle smoothing lengths. 

Similarly, but in spherical geometry, QBisbas et al., 2009] used the HEALPix tessellation ( |G6rski et al., 2005) ) to 
define rays, starting with the lowest level and refining rays into four subrays. Rays are refined when their separation, 
given by the radius r^ay at which they are defined multiplied by the separation angle 9i of the HEALPix level I to 
which they belong, exceeds the local smoothing smoothing length h multiplied by a parameter /2 of order unity. 
Values of /2 of 1.0-1.3 were found to give a reasonable compromise between speed and accuracy. 

Once rays are defined, the discrete integral in Equation[8]is computed along them by defining a series of evaluation 
points, each being a distance fih from the previous one, with fi a dimensionless factor (given a value of 0.25) and 
h being the local smoothing length. A schematic is shown in Figure ID The ionisation front is linearly smoothed 
over one smoothing length and the gas heated accordingly. A smiler adaptive ray-tracing scheme was presented by 
QAbel and Wandelt, 2002| for use on Cartesian grids. This scheme differs from that of [Bisbas et al., 2009] in that 
child rays can be merged in regions where high resolution is not necessary, and that they solve a time-dependent 
problem, using the results of the ray-trace to compute fluxes at cells. The algorithm is taken even further by 
QWise and Abel, 201 1| , who also implement non-ionising radiation (e.g. Lyman-Werner dissociation) and radiation 
pressure. 

QKrumholz et al., 2007b|| use a variant of the ray-tracing method of ||Abel and Wandelt, 2002) , periodically 
rotating the rays with respect to the Cartesian grid to avoid geometrical artefacts. To avoid spurious overcooling at the 
ionisation front, molecular heating and cooling processes are disabled for cells with ionisation fractions in the range 
[0.01,0.99]. They also explore the convergence of the results with changes in the update timestep for the radiation 
scheme, which effectively sets by how much the temperature of a cell near the ionisation front is allowed to change 
in one timestep. Allowing the temperature to change by a factor of 100 led to larger errors in the location of the 
ionisation front at early times, although the error declines as the front expands, and they caution against allowing 
sudden temperature jumps in photoionisation algorithms. This was also pointed out by [Whalen and Norman, 2008], 
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Figure 2: Schematic representation of the adaptive SPH ray-tracing technique employed by QBisbas et al., 2009) . The 
ionising source is represented by the star at the extreme left. Solid lines are rays, with black circles representing the 
evaluation points used to compute the discrete integral Equation!^ Grey circles are locations where they rays are split 
into four sub rays. The dashed line beyond which the particle density becomes abruptly larger is the ionisation front. 


who explicitly compared an algorithm employing the assumption of ionisation equilibrium embodied by Equation [8] 
with a more sophisticated radiative transfer algorithm described in [Whalen and Norman, 20061 . They found that the 
structure of ionisation-front instabilities varied substantially between the two codes, especially at early times, and 
attributed the differences to the sudden heating inherent in the equilibrium method. 


[Peters et al., 2010] implemented the ray-tracing algorithm of | Rijkhorst et al., 20061 in the FLASH AMR 
code ( [Fryxell et al., 2000[ ). The algorithm computes column densities on nested grids using a hybrid long- and 
short-characteristics method. A long characteristic is a ray drawn between a radiation source and an arbitrary cell, 
and may have many segments since it may pass through many intervening grid cells. A single long characteristic 
can transport radiation between two arbitrary points in the simulation. A short characteristic passes only across 
a single grid cell and only transports radiation from one cell to another (see Figure [3 for an illustration). From a 
computational perspective, long characteristics are more amenable to parallel computation, since each ray can be 
treated independently and the radiation transport equation solved along it. However, time is wasted near the source, 
since many rays pass though the same volume. Short characteristics cover the domain uniformly but radiation 
properties of cells must be updated from the source working outwards because each short characteristic must begin 
from the (usuallly interpolated) end solution of one closer to the source. Short characteristic methods are therefore 
difficult to parallelise and more diffusive. 

For a given AMR block, the [Rijkhorst et ah, 2006[ algorithm computes pseudo-short-characteristic rays 
which enter the block from the direction of the source and terminate at all the cell centres for use local to that 
block, and pseudo-long-characteristic ray segments which terminate at the cell corners. It is the latter which are 
shared with other processors so that what are effectively long characteristics can be stitched together across blocks. 
Once the rays have been defined, transfer of ionising photons is solved in a manner similar to that employed by 
[Tremblin et al., 2012b| , save that collisional ionisations are also accounted for. 

Several authors (e.g. [Mackey and Lim, 20101, [Arthur et al., 201 11 ) use the C^-RAY algorithm of 
[Mellema et al., 2006b[ to model photoionisation feedback. This method is photon-conserving and accurately 
solves the RT problem in the (common) case that the computational resolution elements are optically thick. For 
an infinitesimally thin spherical shell of radius r with a radiation source of luminosity L at its centre, the rate of 
absorption of photons per unit area is 
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Figure 3: Illustration of the difference between long and short characteristics. Long characteristics (left panel) are 
drawn from the source to every grid vertex, so that regions close to the source are crossed by many rays. Short 
characteristics (right panel) are drawn across single cells, connecting the point where a vector from the source enters 
the cell (dotted line) to the nearest grid vertex. The radiation held at the entry points of cells (blue circles) must be 
interpolated from values at neighbouring vertices, making this method more diffusive. 


riocal(^ ) 


1 /■°° L^a^exp[T„{r)] 

47rr2 hiy 


(13) 


If, however, the shell has a hnite thickness Ar, photons strike the inner surface at a rate N{r — Ar/2) and emerge 
from the outer surface at a rate N{r + Ar/2), and the number of atoms in the shell is nl^heii- Equation[T3]can then 
be rewritten as the photoionisation rate across a hnite shell. 


1 L^a,^exp[T„{r)] 1 - exp(Arj.)^^ 

47rr2 hv nl4heii 


(14) 


If Ati, tends to zero. Equations [14] reduces to Equation [13] This issue could place severe constraints on the spatial 
discretisation required to model the propagation of an ionisation front. [Mellema et ah, 2006b| also discuss issues of 
time discretisation. The above assumes that the optical depth does not change over the course of the computational 
timestep, so that either a very small timestep is required, or a time-averaged value of the optical depth should be 
used. They show that using at any location the time- and space-averaged values of the neutral density, ionisation 
fraction and the optical depth from the source allows accurate solution using large timesteps, at least as long as the 
local recombination time, and volume elements with large optical depths. 

[Clark et ah, 2012^ present TREECOL which solves Equation ]7] in an SPH code using the gravity tree to speed 
up the calculation of optical depths. The idea behind the gravity tree is that when computing the gravitational forces 
acting on a given particle, groups of particles which are sufficiently far away can be amalgamated into pseudoparticles. 
The gravity tree groups all particles into a hierarchy, usually by recursively dividing the simulation domain into 
eight subdomains. When computing gravitational forces, the angle subtended at the particle by all the tree nodes is 
computed and compared to a parameter 9^, the tree-opening angle. If the node subtends an angle larger than 9c at 
the particle in question, the node is decomposed into its children and they are tested. Once nodes subtending angles 
smaller than 9c are found, they are treated as pseudoparticles. 

TREECOL uses this formalism to save time in computing column-densities along rays. Eor every particle, a 
low-level (48- or 192-wedge) HEALPIX tessellation is constructed and the contribution of tree nodes to the column 
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density in each wedge is computed, respecting the tree-opening angle criterion. This results, for every particle, in 
a moderate-resolution map of the column density from its location to the edge of the simulation domain. This is 
ideal for computing external heating by a radiation bath such as the ambient UV background that permeates the 
ISM. QBate and Keto, 20T5l recently presented a code which hybridises the FLD model of flWhitehouse et ah, 2005| 
with a version of TREECOL to model, respectively, protostellar heating and heating from the intercloud radiation 
field. A method analogous to TREECOL but designed specifically for AMR-based trees was recently described by 
[Valdivia and Hennebelle, 2014| . 


4.1.2 Moment methods 


Many alternative approaches to the RT problem involve so-called moment methods. The fundamental radiative quan¬ 
tity is the intensity or spectral irradiance, I^, which describes at a given location the rate at which energy is emitted 
per unit area, per steradian and per frequency interval. Integrating over frequency and integrating out the angular 
dependence yields the zeroth, first and second moments of the radiation field, better known as the energy density E, 
radiative flux F and the radiation pressure tensor P; 


E = - [ [ i^dndE 

C Ji.=o J 

(15) 

Fi= / Ii/fi.XidHdj/ 

Ju=0 J 

(16) 

Pij = / 4(n.Xi)(n.Xj)dHdi3. 

Ju=0 J 

(17) 


(18) 


While the meaning of E and F are clear, the radiation pressure tensor needs a little explanation. Its (i^j) compo¬ 
nent Pij is the rate at which momentum in the i direction is being advected by the radiation field through a surface 
whose normal is the j direction. 

Once these transformations are done, the radiation field can be treated like a fluid, coupled to the matter den¬ 
sity field by the equations of radiation hydrodynamics (RHD). In a frame following the flow of matter and under the 
assumption of local thermodynamic equilibrium, the RHD equations can be written ([Turner and Stone, 200 1|): 
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(19) 

( 20 ) 
( 21 ) 
( 22 ) 
(23) 


where xf is the frequency-integrated mean opacity including components due to absorption and scattering, and Kp 
and ke are the Planck mean and energy mean absorption opacities, and the colon operator represents a double dot 
product operation. 

A popular approach to solving these equations, owing to its conceptual simplicity, is the flux-limited-diffusion 
(FLD) method. FLD simplifies the evolution of the radiation flux by first assuming a steady state, so that the derivative 
on the LHS of Equation 23 vanishes, then asserting that the radiation held is approximately locally isotropic, so that 

P = E/3 and 

F = -^VE (24) 

3x 
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As the opacity becomes small, this approximation for the radiation flux tends to infinity, whereas in fact F cannot 
exceed cE. In the optically thin limit, the radiation held can in principle be strongly anisotropic, so that the 
assumptions behind the derivation of the moment equations break down. However, this problem can overcome by 
inserting an additional parameter into Eauationl24l 


F = 


3x 


(25) 


where X{E) is iht flux-limiter, whose purpose is to prevent the energy flux becoming unphysically large. The flux 
limiter can be defined via the radiation pressure tensor as follows. P can be written in terms of i? as P = fE, with f 
being the Eddington tensor, which simply encodes the local directionality of the radiation held. Eormally, 


f = ^(1-/)I+^(3/-l)nn 


(26) 


with n = Vi?/1 Vi? I, fin is a tensor formed by the vector outer product of ri with itself, and / is the dimensionless 
Eddington factor. The flux limiter is then defined via / = A + A^i?^, where R = | Vi?|/(xi?)- The function A must 
then be chosen so that the radiation held is always reconstructed smoothly. 

QWhitehouse et al., 2005|| and [Bate, 2009| implemented the effects of accretion heating from low-mass protostars 
in SPH using the ELD approximation. These models account for conversion of gravitational potential energy to heat 
in the accretion flows (i.e. it is radiated away in a physical manner as opposed to be being entirely lost, as in an 
isothermal model, or entirely retained as in an adiabatic model). However, the protostars have no intrinsic luminosity 
of their own, so that the feedback in these calculations is effectively a lower limit ( QOffner et al., 200^ ). 

QKrumholz et al., 2007a| report on a ELD method implemented in the ORION AMR code to model radiative feed¬ 
back in massive molecular cores, considering sink particles as radiation sources. The detailed protostellar models 
are derived from QMcKee and Tan, 2003| and account for dissociation and ionisation of infalling material, deuterium 
burning, core deuterium exhaustion, the onset of convection and hydrogen burning. They set an opacity floor at 
high temperatures, since ELD schemes do not deal well with sharp opacity gradients where the radiation held can be 
strongly anisotropic. 

To avoid the issues that can be encountered in ELD with steep opacity gradients, | Kuiper et ah, 2010) present a 
novel hybrid method which combines ELD and ray-tracing. They split the radiation held into two components. Dif¬ 
fuse thermal dust emission is computed using an ELD method, whereas the direct stellar radiation held is handled by 
doing ray-tracing on a spherical grid either in the grey approximation, or using (typically 60) frequency bins to 
capture frequency-dependent opacities. 

QKrumholz et al., 2009| used the ORION code with ELD to model accretion onto high-mass protostars. They 
simulated a IOOMq, O.lpc radius rotating core which collapsed to form a disc with a central massive object. Once 
the star achieved sufficient mass, Kelvin-Helmholtz contraction raised its luminosity to the point where radiation 
pressure became dynamically important. Radiation dominated bubbles inflated along the rotation axis and infalling 
material landed on the bubbles, travelled around their surfaces and was deposited in the accretion disc. The ac¬ 
cretion rate onto the massive star was thus little altered. A second star grew in the disc resulting in a massive bi¬ 
nary, and the radiation-inflated bubbles became Rayleigh-Taylor unstable, rapidly achieving a steady turbulent state. 

I Kuiper et al., 2012a | and | Kuiper et al., 2012b) simulated essentially the same problem - accretion onto a high- 
mass protostar - using their hybrid ELD/ray-tracing approach and arrived at qualitatively different results. In the latter 
paper, they performed a comparison in which they operated their code using the ELD solver only. Both radiation 
transport schemes drove radiation-dominated cavities, but those produced by the hybrid scheme continued to grow 
until leaving the simulation domain, whereas those from the pure-ELD run collapsed along the rotation axis. They 
found that the cavity in the ELD case was unable to resist accretion onto it, which they attribute to the radiative flux 
in the ELD method tending to point in directions that minimise the optical depth, allowing photons to escape and 
depressurising the cavities. The hybrid scheme does not suffer from this problem, because the stellar radiation held is 
transported by direct ray tracing. They did not observe the kind of instability seen by ||Krumholz et al., 2009) . 

A second alternative to ELD is to compute the Eddington tensor directly. These are usually known as variable Ed¬ 
dington tensor (VET) techniques, and differ in the ways they compute or estimate the tensor. [Gnedin and Abel, 200 1[ 
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( 27 ) 


give a very clear description of the basis of these techniques. They define the VET as /Tr(Py) with 
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^ ij — 



|x-xi|2 


X* — yi\ x-^ — x{ 
|x - xi| |x - xi| ’ 


where xi are the position vectors of the radiation sources and could consist of a set of ^-functions modelling stellar 
sources, or could include every point in the domain if diffuse radiation fields are of interest. T(xi) is the luminosity 
function of the sources and, x is the location at which the radiation field is to be computed. This integral is very 
expensive to evaluate because of the optical depth term 

l IGnedin and Abel, 200 1| make the integral tractable by dropping this term. This means that the radiation flux is 
computed at each point in the simulation domain assuming that the radiation reaching it from all sources suffers only 
geometric dilution, with no opacity effects. This is referred to as the optically thin variable Eddington tensor, 
or OTVET, scheme. They stress that this latter approximation is not the same as in standard diffusion methods, 
because in these the Eddington tensor is computed strictly locally. They also point out that radiation need not be 
propagated at the speed of light, as long as its characteristic velocity is much larger than any dynamical velocities 
present. They suggest that even 100 km s“^ may be adequate. Schemes making use of this approximation are 
known as reduced-speed-of-light (RSL) methods. OTVET methods have also been implemented in SPH codes, 

in gadget-3, and USales et al., 2014) implement an improved version of the 
le in the AREPO code. 

Some tension has recently emerged between ELD and OTVET schemes. [Davis et ah, 2012| use a short- 
characteristics VET method in ATHENA and perform explicitly comparisons with an ELD solver and a Monte Carlo 
solver implemented in the same code. The Monte Carlo algorithm is much too slow to be used in a dynamical 
simulation, so they instead use the three schemes to obtain an equilibrium solution on a single snapshot from a 
shearing-box simulation. They find that the VET and MC methods agree well, with the ELD solver being the odd one 
out, with the discrepancies largest in optically-thin regions. 

[Krumholz and Thompson, 2012) use a two-temperature ELD approximation on a 2D Cartesian grid to study the 
evolution of radiation pressure driven winds in a gravitationally-stratified atmosphere (intended as an approximation 
to ULlRGs and bright, dense young star clusters). Simulations are characterised in terms of whether the radiation 
pressure forces are greater or smaller than the gravitational forces. Where the radiation pressure forces are smaller, the 
gas undergoes vertical oscillations which eventually die out. Otherwise, an instability resembling the Rayleigh Taylor 
instability develops, with columns of dense gas falling into the low-density material at the base of the atmosphere 
where the radiation pressure forces are greatest. These columns contain most of the mass, but the low density and low 
optical depth gas occupies most of the volume, allowing radiation to escape in the vertical direction. The simulation 
reaches a steady turbulent state with nearly constant velocity dispersion and density scale height. 

[Davis et al., 2014| also modelled radiation pressure feedback in ULlRGs. RT was implemented in 2D in the 
ATHENA code using either ELD or VET, and substantial differences exist between the results. In the low-flux 
case, both radiative transfer schemes (and [Krumholz and Thompson, 2012| ) agree that the gas undergoes sta¬ 
ble vertical oscillations. Their high-flux ELD case rapidly becomes Rayleigh-Taylor unstable, as does that of 
[ Krumholz and Thompson, 2012) , and most of the gas sinks back towards the z=0 plane, where it remains in a 
turbulent state. However, in the VET calculation, the behaviour is very different. The RTI also develops, but most 
of the dense gas is nevertheless accelerated upwards out of the simulation domain. The volume-averaged Eddington 
factor in the VET run is generally larger than in the ELD run and exceeds unity for most of the time, while in the 
ELD run it is mostly just under unity. The difference is modest, but crucial. Deeper analysis shows that the two 
schemes agree well in the dense gas, where ELD should be a good approximation, but disagree on the magnitude 
and direction of radiation fluxes in low-optical depth regions where the diffusion approximation is likely to fail. 
In some regions, the ELD fluxes point in the opposite direction to the VET fluxes, accelerating the gas downwards 
instead of upwards, reinforcing the development of low-density channels and accelerating the development of the RTI. 


e.g. I Petkova and Springel, 2009 


I Petkova and Springel, 2009| schei 


4.1.3 Monte Carlo methods 

Monte Carlo methods solve the RT problem by emitting ‘photon packets’ in randomly-chosen directions from the 
radiation source and following their paths through the simulation domain. In each fluid element through which the 
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packet passes the probability of it being absorbed or scattered is computed and more random numbers are used to 
decide its fate. This is repeated with a large enough number of packets to sample the radiation field and iterated 
until some convergence criterion is satisfied. Monte Carlo codes are good at modelling processes such as scattering 
and re-emission, which are cumbersome to compute in ray-tracing schemes, and covering the frequency domain 
of the radiation transport problem is relatively simple, although it usually entails emitting more photon packets. 
The drawback of Monte Carlo methods is that they converge slowly and there is inevitably noise in the resulting 
temperature field owing to the discrete emission of energy. 

Monte Carlo methods are expensive and have traditionally been used to post-process the interaction of non¬ 
dynamic ally-important radiation fields with fixed pre-computed matter density distributions. However, improvements 
in computer power and algorithms have recently led owners of Monte Carlo codes to implement hydrodynamical 
algorithms in the codes (the reverse of what is usually done). 

The SPHRAY algorithm is presented by [ Altay et al., 2008| . Not a Monte Carlo code in the strict sense, SPHRAY 
uses the C^-RAY method to solve the radiation transport equation on randomly-cast rays in an SPH simulation. 
All particles through which a ray passes are identified using the Axis-Aligned Bounding Box method acting on an 
octal tree. In solving the photoionisation problem, the OTS approximation is made (for both hydrogen and helium 
ionisation). The impact parameters of particles intersected by the ray is computed and the smoothing kernel integrated 
through accordingly to compute the particle’s contribution to the optical depth. Photon packets are then propagated 
along rays and a fraction of their energy (1 — e"’’") is subtracted as they pass through each particle. 

[Pawlik and Schaye, 2008| present the TRAPHIC SPH RT code. Radiation packets are emitted from sources every 
simulation timestep and propagated through the gas until a stopping criterion is satisfied. Each source (which can be a 
star particle or a gas particle) emits photons into an array of cones that covers the sky. Virtual SPH particles are placed 
into any cones which do not contain any real gas particles. Photon packets are distributed amongst the real and virtual 
neighbours of a source, and are then passed on in the radial propagation direction in a cone with the same opening 
angle as the original emission cone. The cones therefore subtend smaller and smaller solid angles at the source as 
one moves further away. Gas particles can receive and retransmit multiple photon packets, and packets coming from 
similar directions are merged to improve efficiency. 

Particles absorb energy from photon packets according to their opacities, removing a fraction (1—exp”’’) of the 
energy from the packet. Alternatively, photon packets can be reemitted, treating the gas particle as a radiation source, 
to model scattering. The absorption and reemission process continues until one of two stopping criteria are reached. 
If a state of radiative equilibrium is desired, the process is continued until all packets have been absorbed or have left 
the simulation domain. Otherwise, photon propagation is stopped when the packets have travelled a distance set by 
the speed of light and the timestep. 

I Nayakshin et al., 2009| present an algorithm for modelling radiation pressure in SPH using a Monte Carlo 
method. They track the trajectories of photon packets as r{t) = Tq + Vphot^, where the propagation speed |vphot | need 
not be the speed of the light. They also reduce the packet energies iJphot and momenta Pphot = f?phot/c continuously, 
rather than discretely, using 


1 dpphot 
^phot df 


— PphotPf^' 


(28) 


Packets are destroyed when their momentum drops below 10 of its initial value. 

QHarries, 201 1] [Haworth and Harries, 2012| [H arries, 20T5l implemented a grid-based finite volume hydrodynami- 
cal scheme into the pre-existing TORUS Monte Carlo code. Photon packets have constant total energy and are initially 
given a frequency chosen randomly using the source emission spectrum. The frequency determines the number of 
photons the packet represents. For a source luminosity L, integration time interval At and number of packets N, the 
energy per packet is just e = LAt/N. The packet propagates in a randomly-selected direction for a randomly-chosen 
pathlength I, at the end of which it is absorbed. A new packet is immediately emitted at that location in the same 
fashion, but with a frequency determined by the emission spectrum of the gas at the absorption point. This continues 
until the packet leaves the grid or its propagation time becomes equal to At. 

As packets travel through the grid, they contribute to the energy density in every cell through which they pass. 
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Each packet contributes AU which is estimated from 


AU = 


AttJu 


Av = 
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cAtV 




(29) 


with V being the cell volume. In the photoionisation calculations presented by |Haworth and Harries, 20121 , this can 
be used immediately to solve the ionisation balance equations ( |Osterbrock and Ferland, 2006| ). For a given species 
X with ionisation states X*, etc., the balance equation reads 




1 r°° ATTj^a^{X*)dL' 
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(30) 


where Oi{X'^) are recombination coefficients, a^iX'^) are absorption cross sections, is the electron number density, 
and i/i is the ionisation threshold frequency for species X* respectively. TORUS continues iterating the radiation 
transfer until an equilibrium temperature is achieved in each cell such that the heating and cooling rates balance. In 
radiation hydrodynamics calculations, the radiation transport problem is solved by iteration first, since the solutions 
for subsequent timesteps are likely to be relatively small perturbations on the initial state, and the radiation and 
hydrodynamics problems are then solved one after the other, with the radiation always being done first. 

QHarries, 2015| | discusses two methods by which radiation pressure forces may be computed in Monte Carlo 
schemes. The simplest is to compute the change in momentum Apphot suffered by a packet in a cell and add that 
impulse to the gas in the cell 


Apgas — Apphot — (fiin ^out); (31) 

where Uin and Uout are the unit vectors of the packet’s trajectory as it enters and leaves the cell respectively. The 
radiation force can then be computed at the end of the iteration from frad = X]^Pgas/(AfF). 

The above method suffers problems in optically thin gas, however, where the number of absorptions and 
scatterings can be small or zero. This can be overcome by computing the radiation pressure directly from the radiation 
flux. Equationl^can be rewritten to give an expression for the radiation intensity 


Ii^dildv = 

leading to a Monte Carlo estimate of the radiation flux 


el 
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1 

AtVdv 




(32) 


(33) 


The force may be computed by converting the above expression into one for momentum flux and multiplying it by an 
appropriate opacity 
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(34) 


A much smoother estimate of the radiation pressure is recovered by this last expression, since all packets passing 
through a cell contribute to the estimate, and not just those that are absorbed or scattered. 

In principle, Monte Carlo methods are very easy to parellelise, since each photon packet can be treated inde¬ 
pendently. On a shared memory machine where all the processors can access the entire computational domain, 
parallelisation is than almost trivial. However, most problems are run on distributed-memory machines using the 
message-passing interface (MPI). The MOCCASIN code ( flErcolano et al., 2003| ) gets around this problem by giving 
a copy of the whole domain to each processor, but this is very memory intensive and it is more usual to decompose 
the domain into subdomains, as is done in TORUS. When a photon packet leaves one sub-domain belonging to one 
processor for another belonging to a second processor, the first processor sends an MPI message containing the details 
of the packet to the second. In TORUS, photon packets are communicated in stacks to cut down on communication 
overhead. 
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4.1.4 Alternative methods 


QStamatellos et al., 2007| devised a novel method for estimating the mean optical depth at the location of an SPH 
particle from the local density, temperature and gravitational potential, reasoning that the optical depth and the 
gravitational potential are both non-local properties determined by each particle’s location within the simulation as a 
whole. Each particle is regarded as being inside a spherically symmetric polytropic pseudo-cloud at an unspecified 
location. For any location within the pseudocloud, the central density and scale length are adjusted to reproduce the 
actual density and potential (neglecting any stellar contribution) at the SPH particle’s location within the real gas 
distribution. The optical depth for any given location is computed by integrating out along a radius to the edge of 
the pseudo-cloud. The target particle is then allowed to move anywhere in the pseudo-cloud and a mass-weighted 
average optical depth over all possible positions is computed. Radiation transport is then conducted in the diffusion 
approximation. 

Another approach, taken for example by [Urban et ah, 20091 , is to treat the detailed radiation transport problem 
as subgrid physics, and to parameterise it in some way. [Urban et al., 20091 used detailed one-dimensional DUSTY 
( [Nenkova et ah, 2000| ) models or analytic one-dimensional approximations to compute the temperature distribution 
near protostars owing to their accretion and intrinsic luminosities, and imprint the temperatures on their 3D SPH 
simulations. 

[MacLachlan et al., 2015| use the MCRT Monte Carlo RT code ( [Wood et al., 2004] ) in a novel way to essentially 
post-process the galactic-scale dynamical simulations of [Bonnell et al., 20131 . Snapshots from the SPH simulation 
are interpolated onto a grid and the ionising radiation held from massive stellar sources is calculated using MCRT. 
The solution is mapped back onto the SPH particles and a list of ionised particles generated. The masses of any of 
these that were later accreted by sinks in the dynamical simulation are then docked from the mass of the relevant sink, 
so that the inhuence of ionisation on the star formation rate may be inferred. 


4.2 Winds 


Main-sequence O-star winds have received less attention than photoionisation in the context of simulations of star 
formation, probably owing to theoretical estimates suggesting that photoionisation is likely to be a more important 
feedback mechanism (e.g. [Matzner, 2002| ). There have been many papers written analysing the evolution of wind 
bubbles in ID (e.g. [Castor et al., 1975| , [Arthur, 2012| , | Silich and Tenorio-Tagle, 2013) ), dealing in detail with the 
microphysics at the interface between the hot shocked wind and the cold ISM. However, relatively few authors have 
addressed this problem in 3D. 

Modelling the interaction of stellar winds with the ISM in SPH is difficult because the total mass of the wind is 
much smaller than the mass of molecular gas with which its interaction is to be studied. SPH is most stable when all 
the particles have the same mass. However, this is very difficult to achieve when attempting to model winds, since the 
wind gas may then be represented by too small a number of particles to be adequately resolved. 

Since winds inject momentum as well as matter and energy, [Dale and Bonnell, 2008| took the view that a lower 
limit to their effects could be established by injecting momentum alone. Treating stars as sources of momentum flux, 
they employed a Monte Carlo method in which wind sources were imagined to emit large numbers of ‘momentum 
packets’ in random directions, which were then absorbed by the first gas particle which they struck. A similar 


technique was recently employed by |Ngoumou et al., 20151, although the Monte Carlo element was avoided by 


distributing momentum into the wedges of a HEALPix grid which was recursively refined to ensure that the width 
of the wide end of the wedges was comparable to the local particle resolution at the location where the wind was 
interacting with the gas, in a similar fashion to the ray-casting technique employed by [Bisbas et al., 2009) . 

I Pelupessy and Portegies Zwart, 2012] use the modular AMUSE code to model embedded clusters. Mass and 
energy from winds and supernovae is injected into the gas, which is treated as adiabatic and is modelled using an SPH 
code. Mass loss rates and mechanical luminosities are computed from [Leitherer et al., 1992[ . Gas particles injected 
to model feedback have the same mass as those used to model the background gas and are injected at rest with respect 
to the injecting star. The gas particle masses used are ~ 10“^ — lO^^M© and the wind mass loss rates vary between 
~ 10“® and ^ lO“®M 0 yr“^ per star. At the highest wind mass loss rates, the simulation timestep is such that 
lOs-lOOs of particles are injected per timestep. For the less powerful winds, there are some discretisation issues, but 


18 



































they do not affect the global results substantially. The internal energy carried by the wind particles is determined by 
the mechanical luminosity multiplied by a feedback efficiency parameter which accounts for unmodelled radiative 
losses. The values of the parameter are taken to 0.01-0.1 (see Section 5 for discussion of the results). 

Winds are somewhat more straightforward to model explicitly in grid-based codes and several authors have 
accomplished it, e.g. [Wtinsch et al., 2008) , [Ntormousi et al., 201 1| and |Rogers and Pittard, 20131. Winds are 
modelled by explicitly injecting gas (and therefore momentum and kinetic energy) from the location of the sources 
with mass fluxes determined as functions of time from stellar evolution models. [ Rogers and Pittard, 2013) for 
example simulate the effects of winds on turbulent GMCs using models appropriate for the three O-stars present (with 
initial masses of 35, 32 and 28 M 0 ), with wind terminal velocities fixed at 2 000 km s“^ during the main sequence, 
dropping to 50 km s“^ when each star enters its Wolf-Rayet phase (accompanied by dramatic increases in mass fluxes). 


4.3 Jets and outflows 

The modelling of jets and outflows in simulations of low-mass star formation has become increasing popular in the last 
decade, although the effort has been almost entirely confined to grid-based codes, even though jets can be adequately 
modelled by the injection of momentum and so are less problematic to model in SPH codes than main-sequence 
winds. 

In their fixed-grid MHD calculations, | |Li and Nakamura, 2006| | and [ [Nakamura and Li, 2007) assume that every 
sink particle injects momentum instantaneously into its immediate surroundings. The injected momentum is taken to 
be /M*P*, with / = 0.5 and P = 100km s“^, and M* being the stellar mass. In the earlier paper, the impulse is 
distributed isotropically over the 26 cells neighbouring the cell containing the sink, but in the later work, the outflows 
have two components. Using the direction of the local magnetic held to define the jet axis, they distribute a fraction p 
of the outflow momentum into the neighbouring cells within 30° of the axis. The remainder is distributed as before in 
a spherical component. [Carroll et ah, 2009) adopt a similar mechanism, in which they inject momentum into regions 
with 5° opening angle ten cells across, but with no spherical component. 

I Cunningham et ah, 201 1) describe in detail the implementation of an algorithm to model jets in the ORION AMR 
code, using sink particles as jet sources which inject momentum continuously as the sinks accrete. The mass injection 
rate of the jet is determined by the sink accretion rate in the absence of outflows Mace by Mjet = /w/(l + /w)Macc- 
Conservation of mass results in a modified accretion rate of 1/(1 + /w)Macc, and the jet velocity is set to a fraction /v 
of the Keplerian velocity at the protostellar surface, the protostellar model being derived from [McKee and Tan, 2003| . 
The total momentum injected by a star of mass M* is then /w/vM^^Uk. 

The momentum is introduced over a range of radii Xw between four and eight grid cells from the source falling off 
as r“^, and is modulated by a function of the polar angle 9 and jet opening angle 9q, ^{9, 9q) given by 





(sin)^6» + 6*^) , 


(35) 


derived from [Matzner and McKee, 1999| . The connection with the simulation is made by inserting source terms into 
the density, momentum and energy equations. 

An analogous model is implemented in FLASH by [Federrath et al., 2014| . Outflows are launched in spherical 
cones with opening angle 9q about the sink particle rotation axes, with 9q taken to be 30°. The outflow mass inserted 
in a timestep At is scaled to the sink particle accretion rate so that Mout = fmMAt, with /m taken to be 0.3. The 
outflow mass is inserted uniformly inside the cones, but the outflow velocities are smoothed both with distance from 
the sink, and in an angular sense so that they drop to zero on the surfaces of the outflow cones, avoiding numerical 
instabilities. The chosen smoothing functions are as follows: 
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The outflow characteristic velocity is set to the Keplerian value appropriate for a star with the mass of the sink and a 
radius of lOR©. In order to produce a highly-collimated jet and a more extended ‘wind’, the velocity profile is further 
modified to 

v{0, 0out) = ^0(0,0out) + l^iO, 0out/6). (37) 

The authors stress that they take particular care to correct the mass and momentum fluxes in the two outflow cones to 
ensure that global momentum is exactly conserved. They also transfer a fraction of the accreted angular momentum 
to the outflow. 


4.4 Supernovae 


Few authors have attempted to model the effects of supernovae on well-resolved individual clouds, probably because 
most GMC-scale simulations do not form any O-stars, or do not progress far enough in time for massive stars that do 
form to reach the ends of their lives. However, there are some notable exceptions. 


I Rogers and Pittard, 20131 extend their study of the impact of stellar winds on a turbulent molecular clump 
through the late and terminal stages of their three embedded O-stars. They follow the WR phase of each star and 
allow them to detonate as supernovae one after another by depositing 10®^ erg and IOMq of ejecta instantaneously. 

iDurier and Dalla Vecchia, 2012| set out to compare the relative efficacy of injecting supernova energy in thermal 
and kinetic form in SPH simulations, distributing energy amongst 32 particles nearest the explosion site. Both 
methods reproduced the Sedov solution well when global timesteps were employed, but poorly when individual 
particle timesteps were used. The two methods did agree in the case of individual particle timesteps if they were 
updated immediately after the energy release, and if the timesteps of neighbouring gas particles were forbidden from 
differing by more than a small factor (they suggest 4). 

Using an SPH code, [Walch and Naab, 2014| investigate the detonation of single explosions in clouds of mass 
IO^Mq and radius 16pc, represented by 10® particles. They take the ejecta mass to be SMq and represent the 
ejecta using particles of the same mass as those from which the cloud is built, resulting in there being 80 particles 
in the supernova remnant initially. The ejecta particles are randomly distributed in a O.lpc sphere around the 
explosion source and given a radial velocity of 3 400 km s“^, to give an explosion energy of 10®^ erg (they also 
perform a simulation in which the ejecta particles are not given outward initial velocities, but instead carry the 10®^ 
erg as thermal energy, finding similar results provided small timesteps are used for the ejecta particles and their 
neighbours). Thermodynamics are handled using a constant heating rate and a cooling rate constructed from the 
table in QPlewa, 1995) and the analytical formula in | Koyama and Inutsuka, 20001. Particle energies are integrated 
using substeps in cases when the cooling time becomes shorter than the particle’s dynamical time. They are able 
to accurately reproduce the Sedov-Taylor phase of the remnant evolution, as well as the transition to the radiative 
pressure-driven snowplough stage. 


4.5 Galactic-scale models 


Simulations at the scale of a galactic spiral arm or disc can in general not model most of the feedback processes 
described above for want of resolution. An illustrative problem, for example, comes with the inclusion of supernovae 
in SPH simulations. On galactic timescales, a supernova is the instantaneous point release of a quantity of energy 
which must then be distributed in the gas near the explosion site. The mass resolution of an SPH simulation places 
a strict lower limit on the amount of material over which the explosion energy can be distributed, which in turn 
sets the temperature of the gas. If the mass is too high, the temperature can be much lower than that expected in 
a supernova remnant and is likely to lie in the thermally-unstable regime of the cooling curve. The energy will 
therefore be quickly radiated away. This issue has traditionally been circumvented by temporarily disabling cooling 
at supernova sites. Other problems arise from the tremendous dynamic ranges that need to be modelled. As discussed 
by I Scannapieco et al., 2006| , poor resolution of the interface between hot diffuse material and cold dense clouds, 
particularly in SPH, artificially raises the density and decreases the cooling time in the hot material. 
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I Springel and Hernquist, 2003) implement a feedback model in the GADGET SPH code designed to circumvent 
these problems. They implicitly assume that the gas consists of two phases - a hot rarehed phase and a cold dense 
phase. Mass exchange between the phases occurs via star formation, the evaporation of cold material to form hot 
material, and the cooling and condensation of hot material to form cool material. 

Star formation occurs on a timescale and a fraction /3 of the stellar mass is immediately recycled into the hot 
phase via SNe, so that 


dp^: 

dt 




(38) 


The heating rate from supernovae, expressed in terms of the specihc internal energy stored in the hot phase, is taken 
to be 


d a ^P* 

^(phUh) = esN^=/3usN^, 


(39) 


with esN=10^®erg Mg^. 

The cold phase is assumed to evaporate into the hot phase by a process similar to thermal conduction. The mass 
evaporated from the cold phase is taken to be proportional to the mass released from SNe, so that 


dt ^ 


(40) 


4 

The constant A is environmentally dependent, with its functional form taken to be oc p~'s and the normalisation 
treated as a parameter. 

The cold phase is assumed to grow by thermal instability; 


dpc 

dt 


dph 

dt 


-A(ph,Wh), 

Uh - Uc 


(41) 


where A is a cooling function, and the thermal instability is only permitted to operate in gas whose density exceeds 
a threshold. [Springel and Hernquist, 2003| showed that this model leads to self-regulated star formation, since star 
formation increases the evaporation rate of the cold clouds, which increases the density and cooling rate in the hot gas, 
thereby increasing the rate of formation of cold gas. A reasonable choice of A and the star formation timescale results 
in a star formation law resembling the SK law. 

However, [ Springel and Hernquist, 200^ point out that the assumed tight coupling between the hot and cold phases 
does not permit them to model star-formation driven galactic winds. These are a vital component of galaxy formation, 
since they further suppress star formation by ejecting, or at least cycling, baryons into diffuse regions where star 
formation does not occur, and they also assist disc formation by expelling low-angular momentum material from 
haloes. They therefore additionally implement a parameterised wind model. The wind mass loss rate is taken to be 
proportional to the star formation rate, = pM* and the wind carries a hxed fraction of the total supernova energy 
output. Gas particles are entrained in the wind in a probabilistic fashion, so that in a timestep At, the probability of 
entrainment is 


Pw = 1 - exp 


77(1 — l3)xAt 


(42) 


where x is the mass fraction in the cold phase. In order to prevent wind particles being trapped inside thick discs, 
their hydrodynamic interaction with other particles is disabled for a period of 50Myr. 

[Scannapieco et ah, 2006| improve upon this model with a more sophisticated treatment of SNe. Star particles 
are assigned two smoothing lengths, enclosing equal masses of the hot and cold gaseous phases only. Supernova 
energy is divided between the hot and cold phases, weighted by the corresponding smoothing kernel. The fraction eu 
assigned to the hot phase is injected as thermal energy. A second fraction is assumed to be radiated away by the 
cold phase and lost. The remainder, 1 — Ch — Cr is injected into the cold phase. Cold particles accumulate supernova 
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energy until their thermodynamic properties are similar to their hot neighbours’, at which point they are ‘promoted’ 
to the hot phase. Ch and Cr are then treated as free parameters and adjusted to obtain the desired ISM properties. 

A novel approach to the inclusion of supernova feedback is presented by [Teyssier et ah, 20 T3) , who introduce a 
new variable eturb which represents the density of non-thermal energy and whose evolution is followed via 


Dcturb A PCturb 

P~Rr~ = - 


m 


tdu 


(43) 


where Einj = /0*r7SNesN is the feedback energy injection term composed of the star formation rate, an efficiency 
factor and the energy per supernova, and tdiss is a dissipation timescale, which they set to lOMyr, comparable to the 
typical GMC lifetime (although they discuss several other possibilities). Feedback is connected to the hydrodynamics 
by dehning eturb = pi^turb/^ modifying the total pressure to include thermal and turbulent terms, so that Ptot = 
^therm + po'turb' ^ similar Scheme has been explored by [ Agertz et ah, 2013) . 

Other authors attempt to model supernovae in a more explicit manner in similar fashion to simulations at GMC 
scales. In FLASH QGatto et ah, 20I5| inject 10®^ erg of supernova energy into a region containing IO^Mq of gas or 
eight grid cells across, whichever is larger. If thermalisation of this energy would result in temperatures above 10®K, 
the energy is inserted in thermal form. However, at higher densities and masses, the overcooling problem would 
be encountered, essentially because the Sedov-Taylor phase of the remnant cannot be resolved. [Gatto et ah, 2015| 
get around this problem by computing, given the actual density in the injection region, the expected bubble radius 
and momentum at the end of the Sedov-Taylor phase and injecting this momentum into the cells instead, while also 
heating them to lO^K. 

I Hopkins et ah, 20TTj and [Hopkins et ah, 2012| point out that, on galactic scales, most of the gas is so dense 
that it cools efficiently - only ~10 percent of the total ISM pressure comes from the hot ISM. In denser regions, 
momentum dominates and radiation pressure, stellar winds and supernovae are all comparable when averaged over 
galactic dynamical timescales, j Ostriker and Shetty, 201 1| make a similar point. [Hopkins et ah, 20111 identify star¬ 
forming clumps around the densest gas particles and compute the stellar bolometric luminosity within the clump using 
STARBURST-99 ( QLeitherer et ah, 1999| ) models and a Kroupa IMF. They assume that the momentum is distributed 
equally amongst all gas particles within a clump and for a given particle j in a clump, the imparted momentum flux is 
then 


Pj = (1 + ^pAr)- 


(44) 


with Lj = (Mgas,j/Afgas,ciump)f^ciump- The first factor in the brackets represents momentum deposited in the gas 
by dust absorption of the optical and UV photons from the massive stars. The dust reradiates in the infrared and 
the second term allows for absorption of this radiation, tir is the optical depth across the clump, equivalent to a 
trapping factor, and pp 1 is a parameter which can be adjusted to allow for other sources of momentum, e.g. jets. 


winds and SNe (rjp > 1), or for photon leakage (rjp < 1). A similar model is presented in [Agertz et ah, 2013) but 


their cosmological-scale simulations do not have the resolution to estimate the infrared optical depths, so these are 
estimated from subgrid models. 

As well as the momentum imparted by winds and SNe, they include the thermal energy released by the associated 
shocks, again tabulated from STARBURST-99 models, and including AGB winds as well as main-sequence and WR 
winds. This energy is deposited over the SPH smoothing kernels of the dense gas particles defining the star-forming 
clump centres. Photoionisation heating is implemented by finding particles within the Stromgren spheres centred on 
the clumps and heating the gas inside to lO'^K (or preventing it from cooling below this temperature). 

I Hopkins et ah, 20TT) and [ Hopkins et ah, 20T^ also allow for the fact that feedback partially or completely 
disrupts GMCs, so that substantial quantities of the IR and UV photons released inside them escape and are only 
absorbed at larger distances. | Hopkins et ah, 20T^ allow the energy to spread over the larger of the local gas 
smoothing length or the local gravitational softening length. 

PCeverino et ah, 2014] discuss a similar scheme for including radiation pressure from ionising photons where 
the gas is assumed to be locally optically thin. The radiation pressure is taken to be one third the radiation energy 
density, Prad = 47r//(3c) and isotropic. The intensity / is computed using STARBURST-99 ( [[Leitherer et ah, 1999| ) 
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as a function of the stellar mass, spread over a reference area A, so that I = rm,/A and P^ad = rm,/(i?^c), 
R being set to half a grid cell size for cells containing stellar mass, and one grid cell size for neighbours of such 
cells. If the gas density in the cell exceeds 300cm“^, the radiation pressure is boosted by a factor of (1 + r), where 
T = nceu/300cm“^ to account approximately for the trapping of infrared radiation in optically-thick gas. Using 
CLOUDY ( [Ferland et al., 1998) ) models, they also account of the change in the local heating and cooling rates 
resulting from irradiation by stellar populations of different ages with different SEDs. 

[Roskar et al., 2014) model radiation fields at galactic scales using an escape probability formalism and 
splitting the radiation field into UV and IR components. The energy absorbed from the UV field is taken 
to be i?uv = £^rad[l — exp(—KyvPdustAx)], where Pi.ad=10®^ erg is the total specific energy re¬ 

leased by a IOMq star. It is then assumed that the absorbed UV radiation is reemitted in the infrared, so that 
PiR = Puv[l ~ exp(—KiRPdustAa;)]. The dust opacity is treated as a free parameter. This energy is then added to 
the supernova feedback energy, and effectively deposited as momentum. 


5 What we have learned from including feedback in simulations 

The previous section concentrated on technical descriptions of algorithms and is intended to be mainly of use to 
researchers who are considering writing their own feedback prescription and wish to get an overview of how it has been 
done before. This section is aimed at a different readership and will concentrate on the science results of simulations 
run using the algorithms and codes described. There will inevitably be a small amount of repetition and overlap 
between these two sections, so some forbearance on the part of reader is requested. 

The first five subsections deal with simulations at GMC scales or below. Figure |4] gives an overview of the mass 
and size scales covered by a selection of these simulations. The last subsection deals with simulations at galactic scales 
and above. 


5.1 Fragmentation, the IMF, and star formation rates and efficiencies 


One of the most urgent questions that simulations involving feedback hope to answer is, what is the effect of stellar 
feedback on the star formation process itself? At the smallest scales, as modelled by ||Krumholz et al., 2009|| , 
IKuiper et al., 2012b| (discussed in Section 4.1.2) and, in the case of primordial star formation, by ||Clark et al., 201 1[ 
and QSmith et al., 201 1| , feedback affects the rate at which individual stars build up mass by interacting with accretion 
flows and circumstellar discs, altering their propensity to fragment. At somewhat larger scales, concentrations of 
dense gas will be disrupted, and heating of the gas will raise the Jeans mass and suppresses fragmentation. However, 
shocks driven by expanding bubbles and outflows can also locally increase the gas density and cooling rate. Triggered 
star formation is a very popular topic in observational astronomy, and triggering of star formation by stellar feedback 
is even more intriguing because it should be directly observable in star-forming regions, and gives rise to the attractive 
idea that star formation may be self-propagating. There is a wealth of observational literature on this topic (see 
QDale et al., 2015| | for a recent survey) and it has recently also started to receive increased attention from modellers. 

There are two popular models of triggering - radiation-driven implosion (or cloud-crushing - 
QSandford et al., 1982[ [Bertoldi, 1989| ) and the collect-and-collapse process ( |Elmegreen and Elmegreen, 1978[ ) - 
both of which are now amenable to simulation. 


5.1.1 Radiation-driven implosion 

The RDI or cloud-crushing regards feedback as an external agent which perturbs a stable or quasi-stable equilibrium 
of some kind. It can thus be very rapid and need not involve large masses of material. Many authors have considered 
the effects of irradiating objects of only a few to a few tens of solar masses. Such objects are commonly observed 
around the borders of HII regions, so their evolution is of obvious interest. They are usually modelled as Bonner-Ebert 
spheres. Unless the radiation is strong enough to ionise the entire BBS, it heats a curved layer, thickest at the point 
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Figure 4: Mass-radius parameter space with most of the simulations discussed in this Section overplotted. Colours 
refer to different papers, given in the left and middle columns of the key. Symbols indicate which feedback mechanisms 
are modelled, given in the right column of the key. Lines of constant freefall time (solid) and constant surface density 
(dashed) are also plotted. 


closest to the radiation source. The photo-heated material flows away, driving a shock back into the cloud, a 
phenomenon known as the rocket effect. It is this which may drive the BES to gravitational collapse. 

QGritschneder et ah, 2009al , [Bisbas et ah, 200^ , QBisbas et ah, 201 f| and later [ Ngoumou et ah, 2015) all 
employed SPH simulations and present similar results. [Bisbas et ah, 201 1| performed a parameter study in which 
they varied the mass of the BES and the ionising flux. Low fluxes lead to slow but efficient star formation. Higher 
fluxes produce faster evolution but less efficient production of stars, with star formation concentrated in a pillar-like 
structure behind the ionisation front. This structure is formed by a shock-focussing phenomenon originating in the 
curved outer surface of the BEE. Eor large fluxes or low BES masses, the cloud is destroyed by photoevaporation 
before any stars form. 

I Mackey and Lim, 2010| illuminated groups of triaxial clumps and found that small groups of sufficiently massive 
and dense clumps produced long-lived elongated structures. Three nearly collinear clumps or three clumps close 
together in a triangular configuration were particularly successful in forming pillars by the smearing of the clumps 
away from the ionising source and the Ailing up of shadowed regions behind the clumps with low-density material. 

Other authors have instead turned to irradiating turbulent clouds or boxes. ||Gritschneder et ah, 2009bl and 
iGritschneder et ah, 2010] created a turbulent box which they illuminated with plane-parallel radiation from the 
negative x-direction. The inhomogenous density held allowed the radiation to penetrate to different depths at different 
locations, producing an irregularly-shaped mass of hot gas which then expanded in the face of the ram-pressure of 
the remaining turbulent cold gas. The subsequent evolution was found to depend strongly on the Mach number of 
the initial turbulence. Low Mach numbers presented little resistance to the HII region, which expanded like a piston, 
producing a rather flat ionisation front. Higher Mach numbers allowed progressively longer and more prominent 
pillar structures to project into the ionised gas. In the tips of several of these objects, collapsing cores and discs were 
found, although the simulations could not be run far enough to follow their evolution. 

In a series of papers, UTremblin et ah, 2012b| , ||Tremblin et ah, 2012a|| and ||Tremblin et ah, 2013| thoroughly ex¬ 
amine the irradiation of perturbed ionisation fronts and turbulent boxes with a view to understanding pillar formation. 
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(a) Monchromatic (b) Polychromatic (c) Polychromatic diffuse 


Figure 5: Comparisons of the end result of the medium-flux radiation-driven implosion simulations from 
[Haworth et al., 2(y\2\ . The left panel shows the result of a monochromatic direct radiation field with the diffuse 
field treated using the OTS approximation. The middle panel shows the result of using a polychromatic radiation 
radiation field and again employing the OTS approximation. The right panel depicts the effect of the polychromatic 
radiation field with the diffuse field calculated self-consistently. 


They begin by imposing on an ionisation front a single sine-shaped perturbation of the same length, but different 
widths perpendicular to the radiation field to study the influence of the curvature induced in the photoevaporative 
shock. Narrower perturbations result in longer pillar structures, with shock curvature concentrating material along the 
pillar axis, while a concave pit is excavated around the base of the pillar. They next reintroduce a flat ionisation front 
but with a spherical overdensity just behind it. Shock curvature around the obstacle causes material to meet behind it, 
forming a pillar structure which was observed to be longer the higher the density contrast. Unlike in the case of the 
perturbed fronts, the pillars produced by this process develop pronounced heads shaped roughly like umbrellas. 

In the second paper, a turbulent velocity held with mean Mach number 1, 2 or 4 is irradiated. In common 
with [Gritschneder et al., 2009b[ , they observe that the higher Mach number turbulence is better able to resist the 
compression by the hot ionised gas, and the formation of pillar structures. The pillars have much more complex 
shapes than those formed in the previous paper. In the high Mach number simulation, they also observe globules of 
dense cold gas isolated inside the ionised gas, which are thrown there by the ram pressure of the turbulence. They 
also observe a characteristic double-peaked structure in the gas column-density PDF, with one peak corresponding to 
the turbulent velocity held, and the second to gas compressed by feedback. In the third paper, they show that just such 
PDFs are observed in the neighbourhood of the Pillars of Creation in M16. 

[Haworth and Harries, 20 12^ , [Haworth et al., 20T^ and [Haworth et al., 20T3l examine the subject of triggered 
star formation in bright-rimmed clouds (BRCs) using the TORUS hybrid AMR/Monte Carlo RT code. Their detailed 
treatment of the RT problem allows them to include several physical mechanisms that have been left out of previous 
studies, such as the effect of the diffuse ionising radiation held, and to produce synthetic observational images 
on-the-fly, rather than through post-processing. They And that the diffuse held alters the character of the RDI, 
particular in cases where the radiation held is of moderate strength, where the diffuse held compresses the target 
clump very effectively from the sides, resulting in a much denser and more bullet-shaped configuration, as shown 
in Figure |5] They relate these results to the BRC classification scheme proposed by | Sugitani et al., 19^ , where 
gently-curved rims are classified as type A, tightly-curved as Type B and cometary as Type C. High and low 
radiation fields both produce Type A BRCs, but the strong lateral compression in medium-flux cases reliably leads 
to Type B or C morphology, as shown in Figure |5] Using standard observational techniques on their synthetic 
datacubes, they And that the dynamical states of BRCs can be reasonably well recovered, although the synthetic 
observations systematically underestimate electron densities due to line of sight contamination. This would lead 
to systematic underestimation of the pressure in the ionised flows, and hence to underestimation of the degree of 
shock compression and of the effectiveness of the RDI in real objects. The interpretation of molecular line profiles 
originating in the cold gas is more complex. They are able to reproduce line profiles similar to those observed. 
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but find that some interpretations of those profiles are likely to be erroneous, in particular the two-component 
envelope-expansion-with-core-collapse model invoked to explain the lack of blue asymmetry (which would imply 
infall), lack of any asymmetry, or even red asymmetry (implying expansion or outfiow) often seen in optically-thick 
self-absorbed lines in BRCs. QHaworth et ak, 2013) rule out the EECC model because the material moving towards 
the observer is ionised and not molecular. 

Moving to larger scales, [Dale et ak, 2007bl considered the infiuence of an O-star on a nearly IO^Mq turbulent 
molecular cloud. They used a globally unbound SPH model GMC whose star formation rate and efficiency were 
expected to be low. By comparison with a control simulation, the Lagrangian nature of SPH allowed them to show 
that some of the stellar objects formed in the feedback simulation did not form in the control simulation, but the 
enhancement in star formation rate and efficiency was modest, at 30^0% over sa 0.3 cloud freefall times. The also 
found that some objects which did form in the control simulation were destroyed, aborted or suffered reduced growth 
due to the ionisation of some of the potentially star-forming material. 

[Dale and Bonnell, 2012| performed a similar calculation, except that they irradiated a bound GMC. they found 
that the impact of feedback was even more underwhelming. Although some of the low-density material on the 
outskirts of the cloud was destroyed, and some intermediate-density gas was driven towards the cloud centre by 
the rocket effect, the dense core of the GMC where most of the star formation was taking place remained largely 
unaffected. The number and total mass of stars were little changed and the stellar mass functions proved to be 
statistically indistinguishable. 


5.1.2 The collect-and-collapse process 

In contrast to the RDI model, the collect-and-collapse process is driven internally by stars that have formed inside a 
given cloud and is a large-scale process taking place on relatively long timescales, since it takes time for a sufficient 
mass of gas to be swept up and to become unstable. The mass of gas required is also generally large enough to form 
many stars. 

[Dale et ak, 2007a| used calculations of an HII region expanding in a uniform medium to test a theoretical model 
of the collect-and-collapse process derived by [Whitworth et ak, 1994| from a perturbation-theory analysis of the 
gravitational stability of a shocked gas shell. They obtained reasonably good agreement with the model in terms of 
the time- and length-scales at which the shell became unstable, and of the mass of the fragments produced, and also 
showed that the results were immune to noise of a factors of a few in the initial density field. They found fragment 
masses in the range 10-100 M© but were not able to follow the simulations long enough for many fragments to 
collapse. 

A similar approach was taken by [Walch et ak, 20131 who instead controlled the quantity of structure in their 
initial density field by constructing fractal clouds with fractal dimensions in the range 2.2 (highly-substructured) 
to 2.8 (rather smooth). Clouds with small fractal dimensions (corresponding to large-scale structure) resulted in 
semi-coherent shell structures punctured by large holes through which ionised gas was able to vent (the authors 
refer to these calculations as shell-dominated). Large fractal dimensions, which generate small-scale substructure, 
instead resulted in large numbers of pillarlike-objects pointing towards the ionising source, created by dense clumps 
of material shielding or shadowing gas behind them from the ionising photons (these calculations are hence referred 
to as pillar-dominated). The fractal dimension had a concomitant effect on the fragmentation induced, with low 
fractal dimensions leading to a small number of large fragments and large fractal dimensions producing many small 
objects. The subsequent evolution of the clumps masses is governed by a competition between destruction of the 
clumps by photoevaporation, and their acceleration away from the ionising source by the rocket effect. Very strong 
differences in the clump mass functions result, with the mass function slope being -0.18 for a D=2.0 cloud and -0.91 
for a D=2.8 cloud. Regarding the stars that form, those in the low-D clouds tend to acquire high radial velocities from 
the acceleration of the large coherent clumps in these runs by the rocket effect, so they tend to be found ahead of the 
ionisation front. Conversely in the high-D clouds, the stars are usually found at the tips of pillars and are left behind 
in the HII region. 

[Ntormousi et ak, 201 1| examined a similar process, but acting at still larger scales. They modelled the energy 
and momentum injection by winds of two star clusters 500 pc apart in initially uniform or initially turbulent boxes 
containing 8 OOOK atomic gas. Their intention was to study the formation of molecular material, rather than assuming 
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its preexistence, with the transition from atomic to molecular gas handled by piecewise heating/cooling functions. The 
expansion of the ~ 10®K wind bubbles drove shocked shells into the warm gas which were simultaneously susceptible 
to the non-linear thin shell, thermal, and Kelvin-Helmholtz instabilities, generating very complex structure. The 
former two instabilities combined to form cold overdense clumps in the shells which accounted for 85% of the total 
mass by the end of the simulations after 7Myr. However, the total numbers of clumps was large, so their individual 
masses were equivalent to moderately massive stars. The calculations were rather similar except that additional 
inhomogenieties in the background density field created by turbulence induced the two shells to fragment at different 
times, by approximately 2Myr. The non-uniform emergence of the NTSI also produced larger shear and amplified the 
KH instability, leading to longer filamentary structures in this calculation. Filamentary structures were also generated 
by shear when the two shells collided, resulting a turbulent layer in which the ISM phases mixed and interacted. 


5.1.3 Collapsing turbulent clouds 

Much recent work has concentrated on implementing feedback of one kind or another as additional physics in simu¬ 
lations of collapsing turbulent clouds. These simulations do not fit neatly into the categories of collect-and-collapse 
or RDI, although both processes may be taking place within them at different times and places. However, most of 
these simulations find that the global effect of feedback on these objects is to reduce the star formation rates and 
efficiencies. 

QBate, 20091 use the FLD scheme of QWhitehouse et al., 2005| to model the influence of accretion feedback in 
a 5 OM 0 turbulent cloud. Heating strongly suppresses the formation of new objects after about one freefall time, 
mostly by preventing disc fragmentation. This reduces the numbers of stars formed by a factor of 4 compared 
to a purely barotropic calculation. The decrease in disc fragmentation also results in fewer dynamical interactions, 
sharply decreasing the numbers of brown dwarfs formed. More importantly, the accretion feedback decouples the 
mean stellar mass from the cloud initial Jeans mass. These calculations were extended by [Bate, 2014) to a SOOMq 
cloud to obtain improved statistics. The IMFs produced are statistically indistinguishable from the Chabrier IMF, and 
this result is robust against variations of factors of 300 in the gas metallicity. 

This problem was approached in the ORION code by [Offner et al., 200^ using the FLD implementation described 
by [Krumholz et al., 2007 a| . Protostellar heating again comes to dominate by about one freefall time in the radiation 
transfer calculation. The regions heated by the protostars are small, of order 0.05 pc in radius, and fragmentation 
in most of the cloud proceeds unaffected. The main effect of the feedback is on accretion onto the protostars and 
on their discs. Warmer discs are able to transfer material onto their central stars at higher rates, so can absorb larger 
quantities of infalling gas without becoming unstable and fragmenting. They point out that radiation emitted from 
the protostellar surface originating from, e.g., deuterium burning or Kelvin-Helmholtz contraction, is an essential 
component of feedback in simulations of this kind, and that simulations such as those by [Bate, 2009) which neglect 
it are likely to underestimate the effects of feedback in fragmentation. 

[Krumholz et al., 2010| use the ORION code to investigate how clouds of the same mass, virial ratio and internal 
structure but with different surface densities respond to radiative feedback from accreting protostars. They find that 
low surface density (0.1 g cm“^) clouds analogous to Taurus fragment into a large number of stars of roughly equal 
masses, whereas clouds with surface densities 1 to 10 g cm“^ fragment very little and most of the mass ends up in 
either a single massive binary or a single massive star. The root cause is that the higher density clouds support higher 
accretion rates and therefore higher accretion luminosities, and are also more optically thick so absorb the energy 
released more efficiently. This has the net effect of raising the Jeans mass over substantial fractions of the cloud 
volume, suppressing further fragmentation. 

[Krumholz et al., 201 1| modelled feedback from protostars in a massive (IO^Mq), dense (Ig cm“^) core using the 
prescription from [Offner et al., 200^ , which includes both the energy released by Kelvin-Helmholtz contraction and 
from deuterium burning. By direct comparison with an isothermal calculation, they found that feedback left the stellar 
mass largely unchanged, but led to a smaller number of stars forming. In fact, the warming of the gas eventually shut 
down fragmentation entirely, while allowing accretion onto already-existing stars to continue, resulting in the peak of 
the mass function moving continuously to higher masses. They attribute this to the high density of their clump, which 
allows regions heated by different protostars to overlap, so that virtually all the gas in the clump becomes warm, and 
none of it is able to fragment. 
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In a subsequent calculation, |Krumholz et aL, 2012| tried to remedy this problem in two ways. Firstly, they 
included outflows in addition to thermal feedback. However, this had relatively little impact, reducing star formation 
rates by around 20 percent, entraining very little material and having no major influence on the cloud structure. 
However, they also compared different initial conditions. The earlier calculations were initialised with smooth 
spherical gas distributions seeded with turbulence. Additionally, they constructed initial conditions from fully 
developed turbulence created in a periodic box with no gravity. This produces a density and velocity held which are 
initially self-consistent, leading to more distributed star formation and a slower overall star formation rate, since the 
cloud does not immediately begin to collapse. They found that this lower star formation rate alleviated the overheating 
issue, leading to a consistent mass function at all times. From this study, they conclude that two apparently only 
weakly connected characteristics - the star formation rate and the shape of the mass function - become intimately 
coupled by the inclusion of feedback, and that the latter cannot be correct if the former is too high. 

QFederrath et al., 2014| model the influence of outflow feedback on SOOMq turbulent magnetised clumps with 
virial parameters a = 0.4. They evolve the clumps for almost two freefall times, with a control run used for 
comparison. The feedback simulation forms twice as many sink particles and its star formation efficiency reaches 
75%, whereas that in the control simulation reaches 100%, and the difference in the SFEs is larger at earlier epochs, 
being over a factor of two at one freefall time. On average, the star formation rate per freefall time was 0.57 and 
0.30 in the control and feedback runs. About 40% of the stellar mass in the feedback run was accreted, ejected and 
re-accreted at least once. The average stellar mass is reduced by a factor of RiiS by the combined effect of outflows 
causing more stars to form, and reducing the accretion rates onto individual objects. 

Most simulations of accretion-driven feedback assume that accretion is smooth or continuous. Building on 
several observational studies suggesting otherwise (e.g. [Hartmann and Kenyon, 1996)), flStamatellos et ah, 201 1| 
and flStamatellos et al., 2012| examine the influence of episodic accretion on the fragmentation of circumstellar 
discs. The effects of the energy released by accretion are computed using the radiation transport method of 
flStamatellos et al., 2007| , but the accretion rate is determined using a more sophisticated method than in other 
studies. They divide the accretion discs into two zones - an outer zone, which they can resolve and in which angular 
momentum transport occurs primarily through gravitationally-driven spiral wave formation, and inner disc which 
they do not resolve but instead parameterise, where angular momentum transport is driven by the magneto-rotational 
instability (MRI). The MRI can only operate, however, if the inner disc becomes hot enough to generate the ionisation 
required to couple it to the protostellar magnetic held. The authors show that the ability of the disc to fragment and 
form secondary stars or brown dwarfs is determined by two factors. The first is the ‘base rate’ of accretion through the 
inner disc when it is not experiencing an outburst. This sets the temperature to which the disc (very quickly) relaxes 
when an accretion outburst is complete. The second factor is the duration of the intervals between outbursts, which 
determines whether the disc has time to become gravitationally unstable between stabilising accretion outbursts. 
If the intervals are sufficiently long and the disc is able to cool to low enough inter-burst temperatures, they And 
that accretion feedback is much less effective in suppressing disc fragmentation than inferred by, for example, 
[[Bate, 20091 . | Vazquez-Semadeni et al., 2010) use the ART code ( flKravtsov et al., 1997| ) with energetic feedback 

from O-stars implemented by the deposition of energy in single grid cells where star particles are located, with energy 
deposition rates adjusted to obtain HII regions of reasonable sizes, temperatures and internal velocity dispersions. 
They simulate the formation of a flattened configuration of molecular gas from two colliding streams. Feedback 
increases the mass of dense material while decreasing the star formation efficiency, by factors up to « 10. Feedback 
is unable to destroy the clouds themselves, nor the atomic streams from which they are forming, but it is able to 
destroy the smaller-scale dense clumps where star formation is actually occurring. The HII regions do produce more 
dense clumps, but these generally disperse and fail to form stars of their own. 

The effects of photionizing feedback on embedded clusters formed in artificially-constructed turbulent clouds 
was examined by [ [Dale et al., 2012al and [ [Dale et al., 2013 a) . They found that star formation rates and efficiencies 
were reduced by factors of up to two by the disruption of the dense filaments of gas feeding the clusters. Feedback 
was largely unable to unbind the clusters themselves though and had little effect on the stellar mass functions. They 
were able to demonstrate, by comparison with control simulations, the triggering of stars in the sense of the formation 
of stars that would not otherwise have been born. However, they found that the triggered objects were spatially and 
dynamically mixed with the spontaneously-formed objects and were therefore very difficult to identify. 

Many simulations of turbulent clouds still rely on simple equations of state or optically thin heating and cooling 
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functions to compute the background gas temperature, when in fact this is likely to be set by the external bath of 
radiation and cosmic rays in which the clouds sit, and by the complex chemistry driven thereby. Several papers make 
use of the TREECOL ( [Clark et al., 2012^ ) algorithm to compute the optical depth from the outside of the cloud to any 
point in its interior, and therefore the heating rate from the external radiation field. [Glover and Clark, 2012b| used the 
combination of TREECOL and the chemical network of [Glover and Mac Low, 2007| and [Glover and Clark, 2012a[ 
to evaluate the importance of molecular cooling on star formation. They showed that in fact cooling from C"*" and dust 
were sufficient to allow star formation to proceed, once sufficiently dense gas is formed. They stress that this does not 
mean that molecular cooling is irrelevant, but that it is not necessary. 

[Clark et al., 2012bl use similar numerics to answer the question of how long an observable molecular cloud takes 
to form. They model colliding atomic flows at either 6.8 or 13.6 km s“^. They find that molecular hydrogen appears 
early on, respectively 10 and 3Myr before the onset of star formation, so that the clouds are in this respect ‘molecular’ 
long before they begin manufacturing stars. However, detectable amounts of CO form much later, about 1-2 Myr 
before star formation commences, supporting the idea that there could exist a population of undetectable molecular 
clouds. [Smith et al., 2014] extended this work to galactic-disc scale simulations using AREPO and showed that 42% 
of the molecular mass in their model spiral galaxy was in CO-dark form. 

[Clark and Glover, 2014| examine the question of whether there is a column-density threshold for star formation, 
as suggested by, e.g. [ Schaye, 2004) and [Lada et al., 2010| (who find a value of 116 Mq pc“^. They find that the 
correlation between column- and volume-density in their model clouds is very poor, so that the latter cannot be 
safely inferred from the former, and the star formation rate in volumetric ally-dense gas is much higher than the 
star formation rate in gas at high column densities. They do infer that there is a minimum mean column below 
which molecular clouds are sterile, but is roughly one order of magnitude lower than the threshold discussed by 
[Lada et al., 2010] . 


5.1.4 Star formation from reinserted gas 

Under certain circumstances, the matter injected into star clusters by the combined winds and supernovae of their 
massive stars can itself become the raw material for a subsequent round of star formation. This is a particularly 
intriguing idea, given that many globular clusters are observed to have multiple main sequences. 

Two-dimensional simulations of super star cluster winds were performed by [Wiinsch et al., 2008[ using the ZEUS 
code. Matter and energy are inserted at the centre of a spherical grid. Above a threshold value for the mass injection 
rate and mechanical luminosity, the cluster winds transition from a smooth outflowing state to one in which the matter 
inside a critical stagnation radius is subject to thermal instability. Clumps of gas inside the stagnation radius cool 
catastrophically and collapse under the thermal pressure of the surrounding wind. Most of the collapsing clumps are 
trapped inside the cluster volume and are obvious candidates for forming a second generation of stars. 


5.2 Gas expulsion and cloud destruction 

A long-standing problem in star formation is explaining why it is such a slow and/or inefficient process. The Galaxy’s 
molecular clouds cannot be forming stars on their freefall timescales, because this would result in a Galactic star for¬ 
mation rate about two orders of magnitude higher than is observed, and would have left the Milky Way devoid of gas 
to form stars out of several Gyr ago. Stellar feedback has long been called upon to solve this problem, by slowing the 
collapse of GMCs, or destroying them before they are able to convert more than a few percent of their gas to stars. 

Analytical work by, e.g., [Matzner, 2002[ indicates that expanding HII regions are likely to be the main source 
of energy on GMC scales, at least until the detonation of the first SNe. [Dale et ah, 2005| simulated the impact 
of the Hll region driven by a single very massive star into a non-turbulent cloud. Gravitational collapse gave the 
cloud a filamentary structure, with the filaments meeting at a common hub at the centre of mass, where the massive 
star was to be found. The dense filaments strongly retarded the growth of the Hll region, partly due to the deposi¬ 
tion of neutral gas into the HII region, causing it to collapse and regrow, or flicker. Also, as in the simulations by 
[Gritschneder et al., 2010) , the ram pressure of the flows resisted the expansion of the ionised gas in many directions. 
Much of the ionised gas escaped from the region near the radiation source through moderately collimated outflows. 
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Although sufficient kinetic energy was nominally deposited by ionisation to unbind the cloud, much of it escaped 
through the outflows and collapse continued largely unimpeded. 

In a series of papers, [Dale and Bonnell, 2011) , [Dale and Bonnell, 2012[ and [Dale et al., 2013b| investigated the 
ability of photoionisation to disrupt instead turbulent GMCs with a range of initial radii, masses and turbulent ve¬ 
locity dispersions. The clouds were allowed to form a small number of O-stars, or subclusters large enough to 
host O-stars. The effects of feedback varied strongly with cloud properties, in particular with the escape velocity. 
GMCs in the Milky Way all have very similar column densities, i.e. Mgmc cx i?GMC' Since the escape velocity 
vesc oc {Mgmc / Rgmg)^^'^ , it follows that wesc oc This is weak scaling, but the escape velocities of Milky 

Way clouds lie in the range 1-10 km s“^, which is significant because the speed of sound inside an HII region at 
around solar metallicity is fixed at 10 km s“^. [Dale et al., 2013bl showed using a very simple model that this leads 
to a very steep relation between the mass of a cloud and the fraction of material that can be unbound by HII regions in 
the period before supernova detonation. 

Other authors have investigated the effects of feedback on different geometries. Ionisation feedback form massive 
stars forming inside a rotating IO^Mq clump is studied by [Peters et al., 2010] . The clump contracts to a disc-like 
structure with massive objects forming at the centre. As the simulation progresses, the formation of smaller objects 
further out in the disc starves the massive objects of gas. This in turn allows HII regions to grow, preferentially out 
of the disc plane, although what accretion continues causes them to fluctuate and flicker. They propose that this be¬ 
haviour may explain the well-known ultra compact HII region problem, and find that all the HII region morphologies 
catalogued by [ Wood and Church well, 198^ appear naturally in their simulations. 

[Colin et al., 2013 1 simulate the formation of molecular clouds and clusters in colliding flows of warm neutral 
gas, and follow the effects of feedback from the cluster stars. The clouds formed by their colliding streams are flat¬ 
tened, turbulent, and have masses of order 1O^M0. Ionising feedback is effective at bringing star formation to a halt 
at star formation efficiencies of around ten percent, in contrast to the simulations of [Dale et ah, 2012b| , who found 
that even the lowest-mass clouds were able to continue forming stars, albeit slowly, under the influence of ionisa¬ 
tion. [Colin et al., 2013| suggest that the reason for the discrepancy may be that their clouds, being flattened, are less 
gravitationally bound than those of [Dale et al., 2012b| , and thus easier for feedback to disperse. The fractal clouds 
modelled by [Walch et al., 20121 are also readily destroyed by photoionisation on timescales of 1-2 Myr, despite a 
very inefficient uptake efficiency of kinetic energy of well under 1 percent. Since these clouds are also spherical, this 
explanation cannot work here. It is possible that the fact that the O-stars in [Dale et al., 2012b| ’s calculations first 
have to destroy the dense filaments and accretion flows in which they born impedes them in disrupting the clouds. 
It is also possible, as suggested by [Tremblin et al., 2012a[ , that the turbulence initially present in the cold gas in 
[Dale et al., 2012bl (but not in those of [Walch et al., 20121 ) plays a role. 

While winds are generally regarded as being subordinate to HII regions, they still inject significant quan¬ 
tities of energy into clouds and, in very dense gas, are more efficient at gas dispersal. SPH simulations by 
[Dale and Bonnell, 2008[ and [Dale et al., 2013c| modelled momentum input from O-star winds on turbulent model 
clouds. As with ionisation, they found that the impact of the winds depended strongly on the escape velocity of the 
clouds, despite there being no obvious limit to the rate at which wind bubbles can expand. Winds were able to slow 
the star formation rate somewhat by disrupting the accretion flows feeding stellar clusters, and in general spatially 
separating the stars from the gas. 

As well as the damage they themselves do to GMCs, winds and ionisation are important because of the 
way in which they set the environments in which the eventual supernovae of the massive stars detonate. 
I Pelupessy and Portegies Zwart, 2012) simulate winds and supernova more self-consistently in an SPH simulation 
by injecting hot gas into embedded clusters modelled as Plummer spheres consisting of 10^ stars mixed with 10® SPH 
gas particles whose mass is set to give various SFEs in the range 0.05-0.5, and with total system masses in the range 
7 OO-8OOOM0. Winds and SNe are introduced adiabatically with feedback efficiencies (i.e. fraction of the emitted stel¬ 
lar energy retained by the gas) of either 0.01 or 0.1. For the higher star formation efficiencies, both feedback models 
are able to efficiently expel gas from the clusters, the difference being that with efficient coupling, this is achieved by 
the winds whereas with weak coupling, the supernovae are required to complete the task. Where the SFE is 0.05, the 
clouds are more resistant to the winds owing to the higher gas masses, but are unable to survive the supernovae. In 
these runs, the expulsion of the gas also promptly disrupts the clusters. 

This problem was tackled in an Eulerian context by | Rogers and Pittard, 20131, who modelled the winds (and even- 
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Figure 6: Time evolution of the stellar wind simulation of | Rogers and Pittard, 2013) . The images show density slices 
though the midplane of the simulations 0.67 Myr (left panel), 2.31 Myr (centre panel) and 3.96 Myr (right panel) after 
the ignition of the wind-driving stars located at the centres of the images. 


tual supernovae) from a trio of massive stars in preformed turbulent clouds. The winds rapidly cleared out the central 
region near the stars, but hot wind gas streamed out of the cloud through low-density channels, entraining little neutral 
material, as shown in Figure |6] Much of the coldest densest material was able to survive the winds for several Myr, 
producing pillar-like structures pointing back towards the O-stars. The supernovae did eventually prove sufficient to 
destroy the cloud. Similar models using the NIRVANA code of multiple stellar winds and SNe evolving in a smooth 
background are discussed by ||Krause et ah, 20131 . The energy input from the supernovae is rather unimportant, with 
that accumulated in the wind bubbles beforehand being much more signihcant. The evolution of the supernovae does 
depend on the state of the wind bubble, with larger wind bubbles giving longer cooling times. 

HWalch and Naab, 2014| take the simulations of flWalch et ah, 2012) further by detonating supernovae inside their 
fractal clouds, both those which have suffered ionising feedback and those which have not. The clouds’ fractal struc¬ 
tures distort the expanding shell and allow leakage of hot gas though low-density channels, advecting large quantities 
of energy away. Particularly in the realistic radiative cooling cases, most of the explosion energy is lost and only a few 
percent is imparted to the surrounding cloud. In models where the cloud has experienced photoionisation prior to the 
supernova, the energy uptake by the cloud is more efficient, but only by a factor of ~ 2, because the lower average 
density inside the remnant encounters delays the onset of radiative cooling somewhat. 

RAMSES simulations by Jlffrig and Hennebelle, 2015) investigate the influence of the SN detonation site (which 
other workers hnd to be of crucial importance in galactic-scale simulations - see Section 5.6.2) on a turbulent IO^^Mq 
cloud. SNe are modelled as purely thermal energy injected in a region four cells across, and detonated deep inside the 
cloud, on the cloud border, or outside the cloud. The effect of exploding a supernova outside the cloud was minimal, 
with only a very small fraction of the explosion momentum being transmitted to the cold dense material. Conversely, 
the internal explosion deposits roughly half of its momentum in the cloud and unbinds about half the cloud’s dense 
gas. 


5.3 Driving turbulence 

GMCs exhibit supersonic turbulent velocity helds with characteristic dispersions typically in the range 1-10 
km s“^ (e.g. IHeyeret ah, 20091). The turbulence provides an additional means of support against gravita¬ 
tional collapse and has been incorporated into the very successful gravo-turbulent model of star formation (e.g. 
|Mac Low and Klessen, 20041 ). However, the cause of these velocity helds is much debated, since supersonic 
turbulence, even in a magnetised cloud, will die away due to energy dissipation in shocks in about one crossing time 
( ||Mac Low et ah, 1998] ). Many authors have suggested that, whatever the original source of the turbulence, it is 
continuously replenished by feedback. 

The combined action of multiple jets has been championed in particular by QMatzner, 2007|, who showed 
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analytically that momentum injection from jets could maintain turbulence in ^ IO^Mq clumps. Other authors 
(e.g. [Mac Low and Klessen, 2004) ) note that the rate of turbulent dissipation is comparable to the expected energy 
injection rate from jets. QNakamura and Li, 2007) conhrmed the expectation from QMatzner, 2007[ that turbulence 
driven by collimated outflows decays more slowly than that resulting from multiple isotropic motions. This is a 
consequence of the collimated jets penetrating to greater depths into the surrounding cloud, driving turbulence on 
large length scales, for which the decay times are longer. They observed that even a modest rate of star formation was 
able to reproduce a turbulent-like velocity power spectrum with a power-law form oc k~^'^ over about one and a 
half decades in wavenumber. The turbulent velocity held reached a steady state with gravitational collapse. 

[Carroll et al., 2009] adopted a similar model with 192 randomly-oriented narrow-angle jets inside a periodic 
box, which they also observed to reach a steady state after a few crossing times. They compared the resulting held 
explicitly to an artificially-driven isotropic held with a power spectrum of Ek oc finding that the jets produced 

a spectrum with a slightly steeper exponent of about -2.75. Again, the power-law only extended over approximately 
one decade in k. In these simulations and those of [Matzner, 2007[ , multiple jets rapidly achieve an equilibrium and 
produce a velocity held resembling turbulence, with a power-law power spectrum over moderate ranges of wave 
numbers. However, in both cases, periodic boundaries are used, and it is not clear what affect this has on the results. 
Periodic boundaries make it difficult for injected energy to escape the grid, and makes a steady state a more likely 
outcome. 

[Hansen et al., 201^ explicitly examine the driving of turbulence in their simulations of radiative and outflow 
feedback in low-mass star-forming regions. The action of outflows reverses the decay in the velocity dispersion 
of the gas but the turbulence driven by the outflows is of a different form to the that with which their clouds were 
seeded. In fully developed isotropic hydrodynamic turbulence, the ratio of the energy density in solenoidal versus 
compressive modes is 2, but shear at the edges of their relatively long outflow cavities preferentially drives solenoidal 
modes, increasing this ratio to values of up to 10. Similar results are reported by [Offner and Arce, 2014[ . 

However, not all simulators have arrived at the conclusions that jets are efficient at driving turbulence. 
IBanerjee et al., 20071 modelled in great detail in FLASH the evolution of a single jet and they concluded that the jet 
was not able to propagate the supersonic motions required to drive turbulence to large enough distances, nor was 
it able to entrain sufficient material, for this to be a viable driving mechanism. Jets are in any case not likely to be 
able to maintain turbulence on GMC scales because their combined Ailing factors are too low. Several groups have 
instead looked at the possibility of driving by expanding HII regions, again inspired by analytic calculations (e.g 
[Matzner, 2002| |Krumholz et ah, 2006) ) suggesting that expanding photoionised bubbles should be able to supply 
energy at a high enough rate to compensate for turbulent decay. 

[Mellema et al., 2006a[ simulate the evolution of single HII regions in pre-existing turbulent clouds and And that 
the velocity dispersion in the neutral gas is not only prevented from decaying but raised to values higher than in the 
original velocity held. Similar results were obtained by [Arthur et ah, 201 f] , who also showed that this result is little 
affected by the presence of magnetic fields of realistic strength 

[Walch et al., 20121 in contrast model HII region expansion in fractal clouds with no initial velocities. They 
subtract the kinetic energy due to the radial expansion of the HIIR and And that the remaining random component, 
which they equate with turbulence, is driven more strongly by photoionisation than by gravitational collapse, although 
the efficiency of energy uptake is very low at Ri 0.05%. 

In their simulations of an irradiated turbulent box, [Gritschneder et al., 2009b| measure compressive, solenoidal, 
and total power spectra and compare the simulations including feedback to a control run in which the initial 
turbulence, which has power-spectrum close to Kolmogorov, is allowed to die away. They And significant driving, 
particularly of compressional modes, in the cold gas, but more so at smaller spatial scales, resulting in a flatter power 
spectrum. 

I Krumholz and Thompson, 2012) and [Davis et al., 2014| both And that radiation pressure illuminating an 
isothermal atmosphere from underneath is able to drive turbulence provided that the radiative flux is large enough to 
overwhelm gravity in the atmosphere. This leads to rapidly development of the Rayleigh-Taylor instability and the 
outbreak of turbulence, as discussed in Section 4.1.2. 

Starting from the photoionisation calculations of [Dale and Bonnell, 2012) , [ Boneberg et al., 2015) used structure 
functions and power spectra to assess whether the multiple bubbles expanding from different locations were able 
to replenish otherwise decaying turbulence. The simulations were initially seeded with Burgers turbulence, but the 
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Figure 7: Comparison of the velocity structure function from simulations of turbulent clouds with (red lines and 
symbols) and without (blue lines and symbols) photoionising feedback, 2.2 Myr after the epoch when feedback is 
enabled in the photoionised simulation. The dashed black line is the structure function expected from Kolmogorov 
turbulence. The deviations at large values of dr are due to global expansion of the unconfined clouds. 


velocity field had transitioned to the Kolmogorov slope by the onset of star formation. They found that, particularly 
in the lower-mass and smaller clouds, the structure function in the control simulations subsequently became much 
flatter, with large quantities of power being lost on intermediate scales. By contrast, in the photoionised runs, structure 
functions closely resembling that expected of developed Kolmogorov turbulence were maintained, or restored, as 
shown in Figure|7] Feedback also regenerated the characteristic ratio of power in compressive to solenoidal modes of 
0.5. 

All the work cited in this section so far has concentrated on the study of turbulence in the cold, neutral 
star-forming gas. QMedina et ak, 2014) demonstrated the development of turbulence in the ionised gas inside their 
Fill regions. A statistically-turbulent velocity field takes approximately 1.5 crossing times of the Fill region to 
establish itself, but the power spectrum slopes are substantially shallower than would be expected for compressible 
or incompressible turbulence. They suggest that this is due to the driving of turbulence on all scales, contrary to the 
classical model where it is driven only on the largest scale. 

Many authors modelling feedback on larger scales report producing velocity fields with power spectra re¬ 
sembling turbulence. However, in many cases, feedback is not the only driver of turbulence and it competes 
with shear in galactic discs, self-gravity and various fluid instabilities, particularly the thermal instability. Some 
simulations find that the influence of feedback on turbulence diagnostics such as the gas density PDF is min¬ 
imal (e.g. flWada and Norman, 2007| ). A prominent model of galactic dynamics simulated by, for example, 
[Ostriker and She tty, 201 1|, QKim et ak, 201 1| and [She tty and Ostriker, 2012|, has as one of its main components a 
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vertical equilibrium in the galactic disc between supernova-driven turbulence and turbulent dissipation, which they 
hnd is achieved on timescales shorter than a galactic rotation period. 


5.4 Effects of different feedback mechanisms combined 


Stellar feedback mechanisms obviously do not act alone but in concert and, even if one is likely to dominate, they are 
nevertheless likely to influence each other to some extent. It is often not obvious how such interactions will proceed 
and more and more authors are now investigating the combined effects of several mechanisms. 

Accretion of material by protostars produces two very different but almost always concurrent feedback 
mechanisms, namely accretion heating and jets, which affect the star’s surroundings in very different ways. 

I Cunningham et ah, 201 1| simulated the evolution of a 300 Mq turbulent core including the effects of stellar jets, 
and using FLD to follow the radiative cooling of the shocked gas and radiative feedback from accretion and from 
the protostars themselves. They observe that the outflows reduce the efficacy of radiative feedback by providing 
low-density channels through which photons can escape. 

In common with [Krumholz et ah, 2012|, QHansen et ah, 201^ simulate the interaction of radiative protostellar 


feedback and outflows, using the same model as | Cunningham et ah, 20111. They simulate a 0.65 pc turbulent box 
containing ISSMq of material. Simulations including radiation and outflows form substantially more stars than 
those with radiative feedback alone because the outflows reduce the accretion rates which contribute much of the 
protostellar luminosity. The dual-feedback simulation in fact results in very similar numbers of stars and total stellar 
mass to a simulation with outflows and a barotropic equation of state. This result is difficult to reconcile with that 
of IKrumholz et ah, 2012| , who found that the addition of winds was not able to strongly reduce the protostellar 
accretion rates or luminosities in their calculation. IKrumholz et ah, 2012[ attribute the ineffectiveness of outflows in 
their simulation to the higher densities of their clouds (1 g cm“^ as opposed to |Hansen et ah, 201^ ’s 0.1 g cm 

Another complicating factor in the evolution of combined HII regions and wind bubbles is the possibility of the 
rapid motion of the driving star through its natal cloud. |Arthur and Hoare, 2006| and [Mackey et ah, 20r5) study 
the evolution of wind bubbles inside HII regions driven by runaway stars. The wind bubbles become asymmetrical 
almost immediately, while the HII regions take much longer to do so, owing to their much lower sound speed. 
[ Mackey et ah, 20r5) ’s wind bubbles All about twenty percent of the HII region volume, appropriate for the relatively 
low-luminosity star modelled, and have relatively little influence on the champagne-flow-like internal dynamics of 
the HII region. [Arthur and Hoare, 2006| use higher wind momentum fluxes and the winds dominate the flows in their 
calculations. [Mackey et ah, 20 FS) And, as suggested by [Rosen et ah, 20141 , that much of the wind energy is radiated 
away via microturbulent mixing driven by Kelvin-Helmholtz instabilities at the contact discontinuity between the 
wind and the HII region. However, they stress that the contact discontinuity is likely modified by numerical diffusion, 
and that they include neither physical diffusion nor magnetic fields, so urge caution in interpreting this result. 

From the point of view of individual massive stars, it is clear that photoionisation and the different kinds of stellar 
wind expelled during different stellar evolutionary stages will interact in complex ways. [Freyer et ah, 2003) and 
I Freyer et al., 2006) investigate the coupled evolution of the HII regions and wind bubbles of (respectively) a 60 and 
a 35 Mq star in a uniform background. In the case of the lower-mass object, the effect of the wind is largely to 
compress the HII region into a thick shell lining the inside of the feedback-blown cavity. By contrast, the stronger 
wind of the more massive star sweeps up a shell inside the HII region which becomes dense enough to fragment, 
casting shadows on the ionisation front. This leads to a spikelike structure in the ionised gas, although this does not 
last very long as the HII region is quickly swept up into a thin shell by the wind. 

[Toala and Arthur, 201 1| use 1- and 2D RHD models to study the interaction of the fast WR wind with the 
previously-ejected and slower YSG or RSG wind around 40 and 60 Mq stars, comparing two different stellar 
evolution models. The results depend strongly on the nature of the mass loss in the slower wind phase, on which 
the two stellar evolution models do not agree. A short RSG phase leads to a thin dense wind shell which becomes 
strongly unstable and breaks up into clumps when hit by the WR wind, whereas a long RSGA^SG phase leads to 
a much thicker and lower-density wind which is stable to interaction with the WR wind. Detailed 3D RAMSES 
simulations investigate the mutual interaction between the HII region, wind bubble and eventual supernova explosion 
of a 15M0 star in smooth ambient media at a range of densities. They observe that the expansion of the wind bubble 
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is particularly sensitive to the evolution of the ionising photon flux and wind luminosity of the star as it ages. In 
the denser media, they find that the wind bubble is prevented from expanding by the pressure in the HII region 
until the star moves onto the horizontal branch, and that in general winds from the 15Mq object considered do not 
inject significant energy into the ISM. They also find, in common with [Walch and Naab, 2014| , that the fact that the 
supernova expands into a rarefied medium increases its effectiveness by delaying the transition to the radiative phase. 

The above already complicated scenarios are made richer still by the likelihood that a massive star will have 
a close binary companion, often another massive star. [ Mackey et al., 2014) present a model of the circumstellar 
material around the red supergiant star Betelgeuse. Since massive stars tend to be found in one anothers’ company, the 
cool wind emanating from a red supergiant is likely to be photoionised by its companions. In the case of a star moving 
supersonically through the ISM, the wind will be terminated by a bow-shock, leading to a layered structure where 
a parabola-shaped outer shock encases an ionised wind separated from the supergiant’s neutral wind by a spherical 
neutral shell confined by the ionised wind. The net effect is to confine a large quantity of material (up to about one 
third of the total wind mass) in a dense shell near the star. This shell becomes important when the star eventually 
explodes, since the blast wave will promptly collide with the shock, potentially radiating a large fraction of its energy. 

From the point of view of the effects feedback on the surrounding ISM, flPale et al., 2014| and 
[Ng oumou et al., 2015| studied, at very different scales, the combined effects of photoionisation and momentum- 
driven winds. [ Dale et al., 2014| essentially repeated the simulations of [Dale et ah, 2012bl and [Dale et ah, 2013b| 
with the inclusion of stellar winds, having already established in [Dale et al., 2013c| that winds acting on their own 
were not a very effective means of dispersing any of their model clouds. The results of their simulations were in 
fact little different from the photoionisation-only simulations, as shown in Figure!^ At very early stages, the winds 
helped to clear away the dense filamentary gas in which the O-stars in these calculations are born, but the structure 
of the cold gas rapidly comes to be dominated by the effects of the expanding HII regions. The structures of the HII 
regions themselves were observed to be different, since the winds expel much of the ionised gas through low-density 
channels, excavating large holes inside the HII regions and leaving the ionised gas as a thin skin lining the inside of 
the shocked bubbles of cold gas in a fashion reminiscent of the simulations of [Freyer et al., 20061. 

[ Ngoumou et al., 2015) instead modelled the effect of including winds in the radiation-driven implosion models 
of QBisbas et al., 201 1| . They also observed that, while there are some visible impacts on the morphology of the 
ionised gas, the cold gas is almost entirely insulated from the winds by the HII region. 

[Peters et al., 2014l ’s FLASH simulations include radiative heating from ionising and non-ionising photons, and 
from jets. They start from the same initial state as [Peters et al., 2010] and the collective outflow driven by the stars, 
being directed along their rotation axes, act perpendicularly to the disc, rapidly exiting the simulation box. The 
leading fronts of the outflow cones are Rayleigh-Taylor unstable, exhibiting elongated fingers, and the mass and 
momentum contained in the outflow are substantially larger than those driven by the HII regions. They find in fact 
that only jets are able to reproduce observed levels of mass and momentum entrainment in outflows in their simulations. 


5.5 Combining feedback and magnetic fields 

The role of magnetic fields in ISM evolution remains unclear, partly due to the difficulty of measuring their strengths 
and geometries. Nevertheless, the typical flux densities present in molecular clouds are known to order-of-magnitude 
accuracy, sufficient for their effects to be meaningfully simulated. This is not the place to discuss their general effects, 
but it is clear that jets and bubbles will interact with magnetic fields, and some authors have begun to explore this 
issue. 

[Krumholz et al., 2007b[ considered the effect of a uniform magnetic field on HII region expansion. They 
model the case where the initial gas density and magnetic field strength are such that the thermal sound speed in 
the promptly-ionised gas is much larger than the Alfven speed va in the undisturbed neutral gas, and va is much 
greater than the thermal sound speed in the neutral gas. At early times, the effect of the magnetic field is minimal, 
but it becomes significant when the magnetic pressure inside the swept-up shell becomes comparable to the thermal 
pressure inside the HII region, which occurs when the shell expansion velocity drops to near va- They define a critical 
radius Vm at which this occurs, give by Vm = (cii/uA)^^^f?s- From this point on, the HII region and the shell are 
deformed from the usual spherical shape. Gas motions along the field lines are unaffected, so that the expansion 
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Figure 8: Comparison of the control (top left), ionisation-only (top right), winds-only (bottom left) and dual-feedback 
(bottom right) Run I calculation from QDale et al., 2014| . Colours represent gas column density and white dots are sink 
particles used to model stars. Screenshots are from an epoch 2.2 Myr after the time at which feedback is enabled. While 
the winds do have some effect on the cloud, that of the HII regions is much more severe. This is even more true of 
other simulations in QDale et al., 2014| . 


36 













in these directions continues in a similar way to that in a normal Hll region. Perpendicular to the held, expansion 
is resisted and the shell is initially supported by magnetic pressure. However, as the expansion proceeds, the held 
lines deform to become perpendicular to the shell surface, so the shell loses magnetic support and becomes thinner. 
Very similar results were obtained by [Wise and Abel, 201 1| in testing their RMHD implementation in ENZO, and by 
[Mackey and Lim, 2011) . 

I Gendelev and Krumholz, 2012| build on the work of [Krumholz et ah, 2007b[ by considering also blister Hll 
regions. The orientation of the magnetic held, not surprisingly, strongly affects the results, in particular the degree to 
which the magnetic held lines are distorted and the distribution of the gas kinetic energy in the swept-up shock. In 
particular, orienting the held parallel to the edge of the cloud results in the strongest compression of held lines, and 
in almost all the shell kinetic energy being stored around the lip of the depression in the neutral gas. In all runs, they 
observe that the gas kinetic energy is smaller than in the case with no magnetic helds, but that the total energy is much 
larger, owing to energy stored in the magnetic helds. They propose that this could be an effective driver of turbulence. 

[Peters et ah, 201 1] present an extension to the simulations of [Peters et ah, 2010| in which they combine Hll 
region expansion and magnetic helds in simulations of a IO^Mq rotating cloud. The cloud is initially threaded with a 
uniform magnetic held parallel to its rotation axis, with a mass-to-Hux ratio fourteen times larger than that required to 
support it, consistent with observed held strengths. As expected, while the held cannot prevent gravitational collapse, 
it slows it, reducing the global star formation rate in the calculations by a factor of order unity. However, the most 
massive star, which forms at the centre of the disc, acquires a larger mass because magnetic braking drains the region 
of angular momentum, enhancing infall into the centre. As in their earlier calculations, the central regions of the 
cloud collapse into a battened rotating disc-like structure, with the rotation and collapse dragging the magnetic held 
into a toroidal conhguration. Once a sufficiently massive star has formed to drive an Hll region, the inhuence of the 
magnetic held is largely to slow the Hll region’s expansion. 

[Arthur et ah, 201 1| in contrast examined the effects of magnetic helds on Hll regions expanding in already- 
turbulent clouds. They simulate a 4pc magnetised turbulent box in which the mean atomic number density is 1000 
cm“^ and the ratio of thermal to magnetic pressure is 0.032. Into the centre of the box, they place either an 09 or a 
BO.5 star, whose ionising photon huxes differ by 2.5 orders of magnitude. In the case of the brighter star, their box 
size is smaller than the critical radius dehned by [Krumholz et ah, 2007b[ at which the HII region expansion should 
start to feel the inhuence of the magnetic held. However, in the case of the fainter star, whether or not this is true 
depends on whether the RMS or mean magnetic held strength is used to compute the representative Alfven speed in 
the neutral gas. They hnd that, in the former case and as expected, the presence of the magnetic held has little effect on 
the gas morphology produced by the Hll region expansion, aside from some smoothing out of small-scale structures. 
In the B-star simulation, they hnd that the magnetic held does become critical during the simulation, but that it still 
fails to have a marked effect on the gas structure. They attribute this result to the disordered state of the held, which 
has no preferred direction except on small scales, and cannot therefore deform the Hll region in a systematic or global 
fashion. In fact, because the pressure in the Hll region exceeds the magnetic pressure for most of the duration of the 
simulations, the Hll region expansion deforms the magnetic held into a more ordered conhguration, with held lines 
approximately parallel to the ionisation fronts. 

I Mackey et ah, 2013) consider the inhuence of magnetic helds on the HII regions driven by runaway O-stars 
(using ((-Oph as an example). As in other studies, HII region expansion and the accumulation of dense material 
perpendicular to the held lines are strongly retarded. The opening angle of the cone trailed by the star is increased by 
the presence of the magnetic held. 

The interaction of supernovae and magnetic helds in turbulent IO'^Mq clouds is investigated by 
jlffrig and Hennebelle, 20151, who hnd that the magnetic held has only a rather weak effect, modestly increas¬ 
ing the efficiency of momentum coupling between the supernova and the cold gas in the case where the explosion 
occurs deep inside the cloud. 

Some of the clearest signs of the action of feedback in real astrophysical systems are the production of unusual 
structure, such as pillars. If these are to be correctly interpreted, the role of magnetic helds in their formation and 
evolution must be understood, j Henney et ah, 2009) conduct RMHD simulations of the formation and evolution of 
pillars using the C^-RAY algorithm. Spherical globules with Gaussian density prohles are embedded in an initially 
uniform magnetic held aligned at a given angle to the x-axis, along which lies a point source of ionising photons. 
With magnetic pressures of order 100 times the initial gas pressure, the effects on the evolution of the system are 
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Figure 9: Projections along the z-axis (left) and y-axis (right) of RMHD pillar simulations from 
[ Mackey and Lim, 201 1| at 100 (top) 250 (middle) and 400 (bottom) kyr. The magnetic held is aligned with the 
z-axis. The images are rendered in energy hux in recombination line radiation, arbitrarily normalised. 


profound and depend strongly on the orientation of the held relative to the direction of the source. When the held is 
perpendicular, the globule evolves to a battened structure, whereas a 45 degree angle results in a very asymmetric 
globule which the authors describe as comma-shaped. Thin-shell instabilities caused by rocket acceleration cause 
the heads of the pillars to fragment into smaller objects. Similar results are obtained by [Mackey and Lim, 2011) (an 
example is shown in|9]l, who also observe that radiation-driven implosion and rocket acceleration both tend to align 
the magnetic held with the radiation held. The morphologies in the simulations can be used indirectly to constrain 
the held strengths in real objects, such as the famous M16 pillars. High held strengths produce sheet and ribbonlike 
structures which should be easily observable, but are not seen in M16. Emission in that region is instead consistent 
with that seen in the simulations where the accretion how from the pillar head is roughly spherically symmetric. 

Regardless of their direct interaction with a given feedback mechanism, magnetic helds will always inhuence 
to some extent the environment in which feedback operates. This issue is explored by [Price and Bate, 2009) , who 
extend the work of [Bate, 20091 by repeating their calculations with magnetic helds of various strengths included. 
They hnd that magnetic helds and thermal accretion feedback (which do not interact directly) are additive, in the 
sense that they act at different scales. The magnetic helds, if they are strong enough, affect the large-scale collapse or 
contraction of the cloud, exerting an inhuence on gas of all densities. Accretion feedback, by contrast, only operates 
on the smallest scales by suppressing fragmentation in the densest material. The average stellar mass formed in their 
calculations is lower than in pure hydrodynamic calculations. 

When considering feedback from jets, however, the interaction with the magnetic held is more direct. 
I Wang et ah, 2010) model the effect of jets on magnetised clouds, where the jet emission is aligned with the local held 
direction. Magnetic helds on their own do not much inhuence the denser material from which their sink particles are 
accreting, but the inclusion of outhows partially disrupts the dense gas and slows accretion substantially. Combining 
the two process increases this effect, since the magnetic helds couple gas over much large distances than local pressure 
forces are able to, reducing total mass accretion rates and those onto individual stars by around an order of magnitude. 

[M yers et ah, 2014| perform simulations of the formation of a star cluster including radiative and outhow 
feedback from protostars, and magnetic helds. Contrary to [Price and Bate, 2009| , they hnd that their mean stellar 
mass increases. They interpret this as being due to the mechanism which actually terminates accretion in many cases 
in the [Price and Bate, 2009[ simulations, which is the ejection of sink particles from the densest gas by N-body 
interactions with other sinks, in a globally contracting cloud. This shuts off accretion onto the sink and hxes its 
hnal mass. [ Myers et al., 2014) do not simulate a globally-infalling system, and most of their sinks form from 
fragmentation. As well as slowing the star formation rate by factors of around two, the magnetic helds also suppress 
fragmentation, decreasing the numbers of sinks and raising the average fragment mass. 
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5.6 Large-scale simulations 

In larger-scale simulations, the star formation process itself usually cannot be directly resolved, certainly in terms of 
individual stars. It is often modelled by assuming that all the gas above some threshold density generates stars with 
some efficiency on some timescale, which is often taken to be the freefall time at the threshold density, i.e. 


P{> Pcrit) 

P* = V-r-, -r 

tff (Pcrit ) 


(45) 


With this prescription, the star formation rate automatically reproduces something close to the Schmidt-Kennicutt 
(SK) law. Some simulations do not impose such a law but try to recover it. The SK law is the result of regulation of 
star formation. What is the regulator is unclear, but clearly little can be said about the role of feedback if the SK law 
is imposed from the outset. However, large-scale simulations have other goals, such as reproducing the multiphase 
ISM and driving turbulence, and reproducing the structure of galaxies, especially spirals. 


5.6.1 Vertical slices or columns through a disc and shearing boxes 


Here, simulators model a vertical slice or column through a galactic disc, often including tidal and Coriolis forces as 
external perturbations. These simulations are generally most concerned with the dynamical equilibrium between the 
vertical gravitational field and the thermal and dynamic pressure on the gas. 

The effects of expanding HII regions on galactic disc scales are explored by [ Koyama and Ostriker, 2009) . They 
produce a turbulent, multiphase ISM, although the velocity dispersion of the turbulence (2-A km s is lower than 
would be expected on these scales (7-10 km s“^). They estimate the Toomre Q parameters of their model discs 
and find that they equilibrate into marginally stable states. The star formation rate surface densities also settle into 
equilbria which are mostly well described by power laws resembling the SK law (although they also find this for their 
purely hydrostatic test case, so caution against over-interpretation of this result). 

[Ostriker and Shetty, 201 1) and [|Kim et al., 201 1| discuss a model of galactic evolution in terms of dynamical 
equilibria, where heating balances cooling, pressure balances gravity and turbulent driving balances dissipation. They 
inject momentum from star-forming events and take SNe to be the main sources of feedback. Radiative heating, which 
is assumed to come from massive stars proportionally to the star formation rate surface density, and from a constant 
galactic radiation field, is introduced via a volumetric heating rate. An equilibrium between the vertical gravitational 
acceleration and feedback-driven turbulence is rapidly achieved, heating and cooling and star formation rates settle 
into equilibrium, and thermal and turbulent pressures settle into nearly linear relations with the star formation rate 
surface density. This result, combined with the vertical dynamical equilibrium straightforwardly implies that the star 
formation rate surface density is proportional to the weight of the ISM. 

Similar results are obtained in models of ULIRGs by [ Shetty and Ostriker, 2012) , who include momentum input 
from supernovae. A dynamical equilibrium establishes itself on a timescale short compared to the orbital time within 
the disc in which the input momentum flux balances the vertical gravitational held. The SNe occur only in the densest 
gas in the midplane but are able to accelerate gas to altitudes of more than lOOpc, driving turbulence and resulting in 
an ISM structure dominated by large shells. 

[Kim et al., 2013b| extend the work of [ Shetty and Ostriker, 2012) and [Kim et ah, 201 1| to three-dimensional 
shearing-box calculations. They observe the formation of cold cloudlets by thermal instability which gravitate toward 
and merge with one another in the disc midplane. The models again globally reach a dynamical equilibrium, which 
includes formation and destruction of gravitationally bound clouds. 

Numerous authors have come to the conclusion that the locations and times at which feedback is inserted into 
simulations are at least as important, if not more so, as the detailed physics implemented. | Slyz et al., 2005) simulate a 
three-dimensional periodic non-driven turbulent box 1.28kpc on a side. Star formation is taken to occur in collapsing 
and rapidly-cooling regions whose density exceeds a threshold. Matter and energy released by winds and SNe 
are injected, smeared out over a characteristic timescale of either the local dynamical time, or lOMyr (whichever 
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is smaller). The manner in which the supernovae are detonated strongly affects the resulting star formation rates. 
With the employed time delays, hot gas hlls a progressively larger and larger fraction of the simulation box, driving 
up the density and star formation rate in the cool gas. However, if supernovae are detonated immediately in star 
forming gas, this does not occur and the star formation rates were observed to be up to two orders of magnitude 
lower. More recently, several other groups have explored the issue of the importance of supernova detonation sites. 
QGatto et ah, 2015| obtain radically different results from detonating SNe randomly in their turbulent ISM, or at peaks 
in the cold gas density (arguably more likely explosion sites). Random driving produces an ISM almost entirely 
filled with very hot (10®-10^K) gas with cold material being confined to clumps with small filling factors. Peak 
driving instead produces almost no very hot gas, although there are substantial volumes filled with gas at lO'^K. The 
cold material has a much larger filling factor and a filamentary topology. QWalch et ah, 20141 use the supernova 
prescription of QGatto et ah, 2015| | and examine two additional means of specifying SN sites, namely a mixture of 
random and peak driving, and clustered random driving where the SNe are corTelated with other but not with the 
dense gas. They find that clustering the SNe has only a minor effect. Peak and mixed driving both leave the galactic 
discs strongly concentrated at their mid planes, again with a hlamentary distribution of cold gas. Random driving 
results in an outflow perpendicular to the disc. 

I Hennebelle and Iffrig, 2014] examine the influence of the supernova detonation scheme on star formation in their 
RAMSES models. Peak-driving reduces the star formation rates by factors of 20-30 owing to the rapid destruction 
of dense star-forming material. However, they also find that the details of the implementation of the supernova 
explosions themselves have a strong influence. Introducing the supernova energy in purely thermal form results in a 
star formation rate twice as high as injecting just five percent in kinetic form. 


5.6.2 Whole-galaxy models of spirals 

Although they allow one to achieve better resolution with given computational resources, vertical-slice or shearing 
box simulations require substantial physical components of spiral galaxies, such as spiral arms, to be parameterised. 
Truly satisfying simulations should of course reproduce these features self-consistently. Besides this, self-regulation 
of star formation in the disc plane is a different (though obviously related) question to the vertical equilibria between 
gravity and pressure discussed in the previous subsection. 

In fact, the formation of realistic spiral galaxies, either individually or in the context of cosmological simu¬ 
lations is a long-running problem in computational astrophysics. The issue is catastrophic angular momentum 
loss, in which overcooling of collapsing baryons leads to the formation of too many low-angular-momentum 
spheroidal galaxies, and too few spirals. Early attempts to solve this problem include QKatz and Gunn, 1991] , 
QKatz, 199^ and UNavarro and White, 1994) , which still produced discs that were too small and centrally concen¬ 
trated. QAbadi et ah, 2003] suggested that the resolution of this issue lay in the regulation of star formation by 
feedback. The most important form of feedback at galactic scales has traditionally been considered to be supernovae. 

Notwithstanding these issues, several authors reached the conclusion that many properties of spiral galaxies could 
be reproduced without the need to invoke feedback and that various instabilities, rotational shear and gravitational 
torques were of comparable or greater importance. flWada and Norman, 200 1| And that gravitational and thermal 
instabilities produce a two-phase turbulent ISM in the absence of feedback. After implementing stellar wind and SN 
feedback from Lagrangian star particles, they observe that the overall morphology of their galactic discs and the gas 
density PDF are little changed. Both models produce velocity power-spectra indicative of turbulence, although the 
model without feedback has a steeper slope of -3, compared to -2 in the feedback model, which they attribute to the 
increased importance of shocks in the feedback calculation. Turbulence in both simulations is also driven by shear 
and self-gravity. flKravtsov, 2003) found that their simulations produce a good approximation to the SK law with 
or without feedback. In common with QWada and Norman, 2001| , they And that the effect on the density PDF on 
galactic scales is small, although more low-density gas is produced by feedback. Similar results were reported by 
QWada and Norman, 2007| , despite their rather high supernova rate of 1.5 x 10“^yr“^kpc“^. QRenaud et al., 2013) , 
who include SN and photoionising feedback in their RAMSES models of a Milky Way-type galaxy, concluded that 
feedback was not very efficient in destroying the clouds because of drift between the clouds and the star particles, 
ensuring that most supernovae do not explode inside clouds. Star formation was regulated more by large-scale 
turbulent support than directly by feedback. 
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Other groups have concluded that supernova feedback acting alone can be important, but only if either the 
supernova rate is increased to larger values than observed, as in the aforementioned studies, of if the individual 
supernova energies are increased. [Shetty and Ostriker, 2008) perform two-dimensional simulations, some with 
an external potential to induce the formation of spiral arms, including the effects of SNe modelled as injection of 
momentum to avoid the overcooling problem. Star formation is restricted to gas above a density threshold and SNe 
are introduced probabilistically each timestep. Models with and without spiral potentials both required very strong 
levels of feedback to influence the evolution of these structures, resulting in positive feedback from the collisions 
of adjacent bubbles and a rough equilibrium between cloud formation and destruction. The overall effect is to 
decrease star formation rates, even if positive feedback occurs locally. In AMR simulations of disc galaxy formation, 
I Agertz et al., 2011) included Type la and II supernova feedback, treating the star formation threshold density, 
efficiency and the supernova energy as free parameters. They found that the disc surface density profiles and relative 
contributions of gas and stars to the surface density depended strongly on the parameterisation of the star formation 
rate, whereas the self-regulation of star formation was controlled by the supernovae energies. They were only able to 
achieve self-regulation for explosion energies of 5 x lO^^ergs, rather higher than canonical values. 

Many simulations are, however, significantly influenced by supernova feedback, both in terms of the rates of 
star formation and the structure of the resulting galaxies. These two issues are increasingly found to be connected - 
the distribution of star formation in time controls the final structure of simulated galaxies. [Robertson et al., 2004[ 
employ the | Springel and Hernquist, 2003] multiphase ISM model to model a set of galaxy haloes. One of the 
ten most massive haloes in their simulation hosts a well-resolved disc galaxy with no strong bulge component, 
exhibiting instead an exponential surface brightness profile. Similar results are reported by [Okamoto et al., 2005[ 
and I Scannapieco et al., 2008) who see that feedback suppresses early collapse, leaving large quantities of hot 
gas to cool at late times, from which an extended disc of young stars eventually form. More powerful feedback 
favoured the disc over the spheroid as expected, although the disc-to-spheroid mass ratio never much exceeded 
unity. [Scannapieco et al., 2012) present the Aquila cosmological code comparison project, in which different 
implementations of feedback are tested against each other in different Eulerian and Lagrangian codes, as well as the 
AREPO hybrid code. All the codes implement SN feedback, but using different schemes, although most inject SN 
energy in thermal form. Galaxy morphology is largely independent of which hydro scheme is employed and depends 
instead on the distribution of star formation in time. Forming realistic discs requires that feedback delay star formation 
so that it continues to late times. However, all the codes formed stellar discs which are too concentrated, implying 
a failure to prevent the accretion of low-angular-momentum material. This in turn is attributed to the inability 
of the feedback mechanisms employed to prevent a strong initial burst of star formation. [Aumer et ak, 2013] use 
gadget-3 simulations to study the formation of spiral galaxies in cosmological simulations. Supernova feedback 
is implemented along with radiation pressure, as in [ Hopkins et al., 20TTj . Feedback, in particular pre-supernova 
feedback, helps in production more realistic galaxies by delaying feedback to later times. 

Galactic winds have a role to play in the angular momentum problem, since they preferentially expel low-angular 
momentum material from haloes. Additionally, real galactic discs have much smaller baryon quantities than 
simulations without feedback predict. IGM enrichment indicates that much of the IGM must have been inside halos 
at some point to be enriched, then expelled again. Galactic super winds with mass loss rates much larger than the SFR 
are generally required. 

[Tasker and Bryan, 2006) perform three-dimensional simulations using the ENZO code, with similar star forma¬ 
tion criteria and a similar technique to smear out the effects of Type II SNe feedback in time as j Slyz et al., 2005) . 
The major effect of feedback is to drive a kpc-scale galactic fountain, cycling gas out of the galactic plane, and 
causing the disc in general to puff up. Feedback also increases the star formation timescales so that it is at least to 
some extent self-regulating, but the star formation rates are still substantially higher than observed in real galaxies. 
[ Dubois and Teyssier, 2008) simulate the formation of stellar feedback-driven galactic winds in RAMSES. They allow 
star formation to proceed in the centre of an isolated rotating halo, forming a disc structure. Haloes with different 
masses, spin parameters, virial velocities and star formation timescales are considered, under the influence of SN 
feedback. In the low-mass (10 ^°Mq) haloes, the combined SNe blow bipolar wind cavities, with compact discs 
forming in haloes with low spin parameters being the most efficient wind drivers. In the high-mass IO^^Mq haloes, 
the feedback-driven outflow is stalled by continuing accretion, leading to an accretion shock around the disc, and 
no galactic winds. They conclude that the infalling gas from the halo is a critical determinant in the formation of 
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the wind. However, even in the low-mass haloes, the mass ejection efficiency is about an order of magnitude lower 
than observations suggest it should be, and they are unable to solve the overcooling problem in their simulations. 
[Brook et al., 2012] use the GASOLINE SPH code to study the formation of galaxies over a wide range of masses. 
Early radiative feedback and supernovae are both implemented. Feedback-driven outflows remove low-angular 
momentum material, resulting in more realistic extended discs and removing dark matter cusps. 

Supernovae are generally regarded as the most important feedback mechanism at galactic scales, but several au¬ 
thors have concluded that they cannot on their own reproduce observed galactic properties and that in fact other forms 
of feedback, such as HII regions and radiation pressure do have important contributions to make. [Crain et al., 2015[ 
make the point that the energy, momentum and mass fluxes into the ISM from star-forming regions and supernovae 
are not known, and nor are the magnitudes of the prompt losses of energy by radiation, and of momentum by 
cancellation. In any case, these losses occur at scales too small to be resolved. 

Many groups have tackled this problem by varying the ‘strength’ of feedback in their simulations and ex¬ 
ploring the consequences for observable properties such as molecular cloud properties or the star formation rate. 
[Dobbs et al., 201 ll , for example, simulate whole galactic discs in SPH with a background potential to model 
spiral arms. They include energy from supernovae split 2:1 between kinetic and thermal components, and with 
various assumed star formation efficiencies in dense gas. The higher star formation efficiencies (and therefore 
levels of feedback) produce populations of GMCs which are predominantly unbound and irregularly shaped, 
although they note that collisions between clouds also contribute substantially to their internal velocity dispersions. 
[Dobbs and Pringle, 201 3j quantify this result further and find that it is mostly the smaller IO^^Mq) clouds that are 
strongly affected by feedback, whereas the larger ^ IO^Mq are more affected by shear. [Dobbs, 2015| resimulate 
portions of | Dobbs and Pringle, 2013 |’s calculations at higher resolution and explore different prescriptions for the 
timing of the injection of supernova feedback, including delayed and stochastic models and smooth and instantaneous 
ene rgy injection, but th ey find that in practice their choices have little impact on the outcome of their simulations. 

i tibler et al., 2014) explicitly test the effect of varying the strength of combined-feedback models. Their strong 
feedback model is derived from [Aumer et al., 20131 and includes SNe kinetic and thermal feedback in a multiphase 
decoupled ISM as described by | Springel and Hernquist, 20031, as well as radiation pressure, while their weak 
feedback model includes only SNe. The star formation efficiency of the weak-feedback model is always too high, 
and the resulting galaxy is spheroidal, whereas the strong feedback scheme gives disc-like stellar distributions and 
credible baryon conversion efficiencies at z=2-3, but too few stars at z=l. Higher star formation rates in the weak 
feedback case also result in older, redder stellar populations. 

[Ceverino et ah, 2014| improve on the basic thermal feedback from winds and SNe in the ART code 
( [Kravtsov et al., 1997| ) by implementing radiation pressure feedback. They find substantial differences, with 
radiation pressure driving additional turbulence in galactic discs and destroying the most compact and high-column 
density clumps, which are similar to normal GMCs and substantially less massive than the giant clumps of the kind 
observed in moderate redshift galaxies (e.g. [Genzel et al., 20lT| ). The stellar distribution is correspondingly more 
extended and less concentrated, and star formation rates drop by factors of a few. The formation of smaller satellite 
galaxies is almost entirely terminated. 

The potential of feedback to reduce the over-formation of stars in cosmological SPH simulations is explored 
by [Stinson et al., 20l'3l , who found that supernova feedback on its own is not sufficient and turn to pre-supernova 
radiative feedback. They do not include kinetic feedback as, for example, in [Hopkins et al., 20TT) , because molecular 
clouds in their simulations are not resolved, being represented by only a few particles. They instead introduce thermal 
feedback over a 4 Myr period into these few particles. They do not disable cooling, so the gas rapidly cools to lO^K. 
The large-scale dynamical effects are negligible, but star formation is locally suppressed. They adjust the efficiency 
of energy uptake to reproduce the required stellar- to halo-mass ratio, and find that the resulting star formation rate is 
highly sensitive to this parameter. 

A suite of AMR simulations of an isolated disc galaxy is executed by [Agertz et al., 2013) in which the feedback 
prescription is varied. Momentum is injected either by momentum kicks or as a non-thermal form of pressure 
and different values of the IR optical depth for momentum absorption are tested, cooling delays of 10 or 40 Myr 
after energy injection are implemented, and a passively-advected decaying energy variable akin to that used by 
[ Teyssier et al., 20T3j is also explored. They conclude that pre-supernova feedback is a crucial ingredient in their 
models, capable of reducing star formation rates and efficiencies by an order of magnitude on its own, and averting 
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the overcooling problem by destroying dense gas before SNe can detonate inside it. This very rich suite of feedback 
models was used by [ Agertz and Kravtsov, 2015] to show that feedback is only effective in regulating star formation 
rates if local star formation efficiencies are high enough that there are strong spatial and temporal correlations between 
the injection of momentum and energy, but that in such cases, the results are sensitive to the details of the feedback 
prescription. 

[ Moody et al., 2014| test two models of feedback in their simulations of the formation of galaxies, and in 
particular of the very massive (10^-10®Mq) and extended (~ Ikpc) clouds observed in star-forming galaxies at 
redshifts z ~ 1 — 3 (e.g. [Genzel et al., 20lT| ). The first model encompasses only thermal feedback from stellar winds 
and SNe, modelled as a constant heating rate for 40Myr after a star formation event. The second model includes 
radiation pressure from the ionising radiation fields of the sources, suppressing the formation of lower-mass clumps, 
and slowing star formation by factors of around two. The addition of radiation feedback produces a much stronger 
decrease in the star formation rate in their simulations than SN feedback on its own, but they still find that they 
overproduce stars by factors of a few. They also find that stellar feedback reduces the bulgeidisc ratio and produces 
discs with flat rotation curves, but which are too thick and only poorly rotationally-supported. 

I Schaye et al., 2015| imple ment star formation stochastically via a scheme proposed in 
I Schaye and Dalla Vecchia, 2008 J in which the rate is set by the local pressure, and which automatically repro¬ 
duces the SK relation. Star formation events inject energy and momentum from radiation, winds and SNe. Energy is 
injected via a temperature increment to avoid spurious overcooling, with an efficiency factor which sets how much 
energy is available for feedback. They adjust parameters such as the density threshold for star formation and the 
efficiency of energy uptake, calibrating the results using the z = 0.1 galactic stellar mass function as a yardstick. 
Since feedback regulates the supply of gas to the ISM, the strength of feedback influences a wide variety of properties 
of the formed galaxy population, including central black hole masses, spatial concentration of stars, star formation 
rates and metallicity. 

Other authors, in particular [Hopkins et al., 20TT| , [Hopkins et ah, 20T2\ , [Hopkins et al., 20T3) and 
[ Hopkins et al., 2014) attempt to model feedback in a physically-motivated manner with as few adjustable pa¬ 
rameters as possible. These papers model isolated and merging galaxies of several different types - high-redshift 
objects, and local Milky Way-type and dwarf galaxies. The structures of the merging galaxies are largely controlled 
by gravitational torques, and the violently unstable high-redshift galaxy is less affected by feedback and more by 
gravitational instability. However, the ISM in their model galaxies generally reaches a steady state after a few 
dynamical times in which the star formation rate and feedback keep the galactic discs on the edge of Toomre 
instability, as observed by [ Koyama and Ostriker, 2009) . Good approximations to the SK law emerge naturally from 
this state. For the larger galaxies, radiation pressure is the dominant feedback mechanism, particularly in driving 
large-scale outflows. However, the density of their dwarf model is low enough, and the corresponding cooling time 
long enough, that thermal feedback from SNe and winds become important. [ Hopkins et al., 20T4) conclude that 
feedback is both necessary and sufficient to regulate star formation, but only when all their feedback mechanisms 
act in concert. The results of these simulations are an interesting contrast to those of [ Ostriker and Shetty, 201 1) 
and [Kim et al., 201 1| , who discuss vertical stability maintained by an equilibrium between gravity and dynamical 
pressure acting perpendicular to the disc plane. 

[Roskar et al., 2014) study disc galaxy formation in RAMSES and examine the importance of pre-supernova 
radiative feedback as well as supernovae. They vary the dust opacity so that the gas absorbs larger and larger quantities 
of momentum form the stellar radiation held. They And that radiation feedback is able to avert the overproduction 
of stars by expelling gas from their model halo. However, they also And that this result cannot be achieved without 
severely disrupting the morphology and dynamics of the disc. 


5.6.3 Dwarf galaxies 

Spiral galaxies have received a lot of attention from modellers, not least because we live inside one. However, 
dwarf galaxies are the building blocks of larger systems, are the most dark-matter dominated objects known, and are 
enigmatic in the sense that cosmological simulations produce many more of them than are actually observed. They 
are thus a rich held of study, and their smaller masses makes the effect of feedback upon them more dramatic. 

Many authors have addressed the self-regulation of star formation in dwarf galaxies which can, under the influence 
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of feedback, operate in the the so-called breathing mode ( [Stinson et al., 20071 ), where gas is expelled into the halo 
but eventually reaccretes to form a new generation of stars, which restart the cycle. This is described in detail by 
[Kawata et al., 2014) , who perform SPH simulations including stellar winds and SNe, with radiative cooling disabled 
for the duration of feedback. Star formation in the feedback simulation becomes intermittent, and large fractions of 
the galaxy volume are cleared of gas. The expanding bubbles driven by the earliest-forming stars collide with gas 
still accreting onto the galaxy and a second round of star formation is induced. Feedback from the second generation 
of stars increases the bubble size further, temporarily terminating star formation. However, the bubbles are not able to 
expel gas entirely from the galaxy’s halo and it eventually reaccretes for yet another round of star formation. 

[Kim et al., 2013a| model dwarfs using the ENZO code, injecting matter and energy from SNe with a characteristic 
time delay of 6Myr, and including transport of ionising radiation using the methods of [Wise and Abel, 201 1[ and 
[Kim et al., 201 1| , with star particles given ionising fluxes proportional to their masses during the 6Myr interval 
before the SNe. Both heating and momentum from the radiation held are included. They observe that star formation 
becomes self-regulating after a transient period of order lOMyr, and that the addition of radiative feedback slows the 
star formation rate by roughly one fifth compared to models involving only SNe. 

Some dwarf galaxy simulations reveal interesting parallels with simulations conducted at much smaller molecular 
cloud scales. [Pawlik et al., 20131 use the TRAPHIC code to investigate the formation of the first galaxies under the 
influence of radiative feedback from the first stars. The primary coolant in the metal-free gas is rovibrationally 
excited molecular hydrogen, which is destroyed by Lyman-Werner band photons, so that stars deprive forming 
galaxies of their ability to cool. Secondly, hydrogen-ionising photons heat the gas, driving it out of haloes with 
virial temperatures lower than lO'^K, and slowing accretion onto lower-mass objects. However, photoionisation also 
produces free electrons, catalysing molecular hydrogen formation and increasing the ability of the gas to cool. They 
And, as in smaller-scale simulations, that the effects of radiation are constrained by inhomogenous gas structures, 
in this case filamentary accretion flows feeding baryons into their simulated forming dwarf galaxy. The similarity 
between this and the process observed by [Dale et al., 2005[ is striking. The reduction in accretion is sufficient to 
terminate star formation, eventually enabling the gas to cool and to undergo another burst of star formation, by 
which time the halo is too massive to lose any further baryons due to feedback, and the second starburst lasts longer. 
The virial temperature of the halo becomes high enough to enable atomic hydrogen cooling and the galaxy forms 
stars continuously, although more slowly than in the case with no radiative feedback, anf it is the photoionisation 
heating and not the Lyman-Werner photodissociation that makes the difference. As the simulation proceeds, the gas 
baryon density in the halo rises and the expansion of the HII region is retarded by the high recombination rates (as in 
[Dale et al., 2005[ and [Peters et ah, 2010] ), at which point the simulation effectively ceases to respond to feedback. 
Radiative feedback also effects, via gas expulsion and redistribution, the dark matter distribution, decreasing the 
central density. The structure of the baryonic discs is largely unchanged, although photoheating suppresses star 
formation in neighbouring haloes as well. 

Other authors model situations in which star formation is overall slowed or outright terminated by feedback. 
[Wise et al., 2012| , who consider radiation pressure from population III and II star formation, implemented using the 
scheme of [Wise and Abel, 201 1[ . They focus on the evolution of the most massive halo in their simulation, finding 
that radiation pressure reduces the overcooling problem and reduces star formation by a factor of Ri 5. Thermal and 
momentum transfer from the stellar radiation held create a ~ lOOpc expanding turbulent region centred on the main 
concentration of stars. 

In many of these models, external feedback is at least as important if not more so than internal processes. 
[Whalen et al., 20081 examine the external irradiation of minihaloes by massive primordial stars. Their results show 
interesting parallels with the much smaller-scale, present-day calculations of, for example, [Bisbas et ah, 201 fl . 
They show that radiative feedback accelerates star formation in relative dense and evolved haloes, but terminates it in 
more diffuse structures. 

[Petkova and Springel, 201 f[ used their optically-thin variable Eddington tensor radiation transport scheme to 
examine the feasibility of cosmic reionisation by star-forming galaxies, and the consequences on the global star 
formation rate. Photo-heating of the intergalactic medium significantly slows the accumulation of baryons in dark 
matter haloes, particularly for the lower-mass haloes, and slows star formation by a factor approaching two. Although 
they are mainly interested in the effects of supernovae, [Sawala et al., 2010) And that background UV radiation is 
of crucial importance in their simulations. The combined mechanisms are sufficient to empty dwarf haloes of gas 
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and terminate star formation by z=6. | Simpson et al., 20T3) reached similar conclusions using ENZO. Their halo’s 
quantity of dense star-forming material is largely regulated by SN feedback, but the total gas fraction is largely set by 
photoheating, and once it is cleared of gas, it cannot reaccrete. 

Reionisation and photoheating is also examined by flPawlik et ah, 20151 , who point out that reionisation is 
self-regulating. Photoheating boils gas out of low-mass haloes and shuts down accretion, lowering the ionising 
emissivities of galaxies. However, a positive feedback is also exerted, since the heating smooths out density 
inhomogenieties and reduces the recombination rate in the IGM. SN feedback also affects reionisation, since it 
clears out gas from low-mass haloes, shutting down star formation, although it also potentially locally triggers star 
formation and opens low-density escape routes for photons, which may enhance reionisation. Star formation is 
suppressed by almost an order of magnitude by feedback between z = 14 and z = 6, with SNe being much more 
significant than photoheating, although the two mechanisms complement each other, suppressing star formation more 
when acting together than alone. SNe suppress star formation evenly across most of the range of galactic masses, 
whereas photoheating suppresses star formation only in the lower mass haloes. In both cases, suppression occurs by 
destruction of dense gas and by reducing accretion. 

The internal structures of dwarf galaxies are also sensitive to feedback. The core-cusp problem in dwarf galaxies 
is discussed by | Teyssier et al., 20131. While pure N-body models predict that dark-matter haloes have central cusps, 
observations suggest that they instead have constant density cores. Of many possible solutions, the least outlandish is 
the modification of the baryon distribution and hence the dark-matter distribution by feedback. They examine this 
question with the RAMSES AMR code. The effect of feedback is to produce a thick, turbulent gas disc supporting 
outflows and emptying the central regions of baryons. The star formation rate is lowered, on average, by an order of 
magnitude and is much more staccato in nature. Gas cools and collapse, starting a burst of star formation which expels 
gas from the disc, shutting star formation down. However, gas falls back onto the disc and triggers another round of 
star formation, and the cycle repeats. The effect of gas removal form the central regions is indeed to transform the 
dark-matter distribution there from a cusp to a core after « 0.5Gyr. 


6 Summary and outlook 

The last decade has seen the development of a tremendous variety of sophisticated algorithms to model the various 
kinds of stellar feedback, and a corresponding wealth of simulations employing these algorithms to answer a wide 
variety of questions over a huge range of scales. We have learned an enormous amount from these simulations. 

Feedback regulates or helps to regulate the rate at which gas is converted to stars at every stage of the star 
formation process. The background cosmic ionising radiation field controls the accumulation of baryons in primordial 
haloes, and supernova and radiation pressure feedback are capable of emptying haloes of gas and terminating star 
formation inside them (e.g. [jSawala et al., 2010| ). Most calculations agree that the cycling of baryonic matter between 
the hot and cold phases of the ISM and the formation of GMCs is influenced - if not controlled by - feedback. 
Other mechanisms, such as gravitational torques, are of course involved and the relative contributions of the various 
processes is still a matter of debate. It seems to be increasingly clear that, whether they are dominant energetically or 
not, SNe are not the only important form of feedback and other mechanisms, particularly radiation pressure, cannot 
be ignored. This is particularly true in simulations of dwarf galaxies, whose lower escape velocities make them 
vulnerable to radiative feedback (e.g. flSawala et al., 2010) , flPawlik et ah, 20l5| ). 

However, the majority of galactic-scale simulations are still not able to resolve these processes. Some authors 
have parameterised the strength and form of feedback and varied the parameters until acceptable fits to some 
observable metric(s) are obtained (e.g. flSchaye et al., 2015| ). A more satisfying approach, taken for example by 
[ Hopkins et al., 20l4) , is to try to devise physically-motivated subgrid models, but even these must rely on some 
physics, such as the leakage of ionising photons, which cannot be resolved in the simulations themselves. 

Simulations at GMC scales have the advantage that they have much better length and mass resolution and have 
shown that all forms of feedback play some role in regulating the rate and efficiency of star formation in these objects, 
from radiation pressure on accretion flows at sub-AU scales to winds, HII regions and supernovae up to ~ lOOpc 
scales. However, most simulations still overproduce stars and none are yet capable of terminating star formation and 
reaching ‘completion’. 


45 



















These models are capable of realistically reproducing the geometrical structure of clouds and therefore also 
quantities which depend on this, again such as photon leakage, can be much more accurately measured. So far, almost 
no effort has been made to connect simulations performed at these smaller scales to galactic-column, galactic-disc 
or cosmological calculations. The GMC-scale models are presently an untapped resource which could supply 
more accurate parameterisations of many quantities of use in the larger-scale simulations. However, simulations by 
QDale et al., 2014| suggest that the permeability of clouds to photons, momentum, energy and polluted ejecta, is a 
cloud-dependent property, making its parameterisation more difficult. 

In addition, none of the GMC-scale simulations yet include all feedback modes. It is often said that HII regions 
are likely to be the most important feedback mechanism on GMC scales, at least until the detonation of the first 
supernova. While this is probably true, it does not mean that other types of feedback can be ignored. Some molecular 
clouds, such as Ophiuchus, are too small to form any OB stars and are of necessity dominated by accretion feedback. 
Such clouds are the most common by number and, in galaxies such as M33, which has a very steep GMC mass 
function, they also dominate the molecular mass. QOffner et al., 200^ makes the point that even in clouds that are 
massive enough to manufacture O-stars, accretion is still the dominant mode before these massive stars are born (and 
may continue to be even afterwards in places which are shielded from ionising photons). Accretion feedback therefore 
does help determine the environments in which the massive stars form, particularly cloud properties such as the star 
formation efficiency at that epoch. This can also be said with regard to magnetic fields. 

Two problems which are therefore of crucial importance are determining how all the different feedback types 
interact with one another in clouds with various different properties (density, escape velocity, geometry, etc.), and 
how magnetic fields contribute to this picture. It is clear from work already done that different feedback mechanisms 
are not necessarily additive (e.g. | Myers et al., 20141), and that the likely most important mechanism - Hll regions 
- can be strongly affected by magnetic fields (e.g QGendelev and Krumholz, 2012| |). In addition, the work of 
I Hennebelle and Iffrig, 2014) shows that the effects of SN feedback on the largest scales is likely to be strongly 
dependent on the details of the environment in which the massive stars explode, which sets the relative quantities of 
thermal and kinetic energy deposited. Teasing out all these interactions requires a great deal of painstaking work, 
particularly if we wish to explain, as opposed to simply reproduce, the evolution of molecular clouds. 

They are some important disagreements between models that need to be resolved. Those that have emerged 
between flux-limited diffusion and variable Eddington tensor radiation transport methods are of particular concern. 
However, the question of why some galactic disc simulations require feedback to produce realistic galaxy properties 
whilst some can achieve them without feedback (e.g [Wada and Norman, 2007[ ) is also curious and needs to be 
explored, as is the differing opinions of, for example, |Colm et al., 2013| and QDale et al., 2014| on how efficiently 
ionisation is able to disperse intermediate-mass GMCs. While much progress has been made, and the pace is 
accelerating, many of the details of the effects of feedback on star formation, the ISM and galactic structure are still 
murky. 
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